
# Libraries

install.packages(c("plyr","psych", "dplyr", "ggplot2", "weights", "mediation", "scales", "jtools", "ivreg", "stargazer", "performance"))

library(plyr)
library(psych)
library(dplyr)
library(ggplot2)
library(weights)
library(mediation)
library(scales)
library(jtools)
library(stargazer)
library(ivreg)
library(performance)


set.seed(89)


#######################
#### Load and rename datasets 
#######################

load("Redistribution_Datasets.RData") #Make sure that the RData file is saved in the same directory as your project to enable easy loading. 

#Renaming country datasets
au1 <- au_wave1
ch1 <- ch_wave1
de1 <- de_wave1
fr1 <- fr_wave1
uk1 <- uk_wave1
us1 <- us_wave1
au3 <- au_wave3
ch3 <- ch_wave3
de3 <- de_wave3
fr3 <- fr_wave3
uk3 <- uk_wave3
us3 <- us_wave3


#######################
#### Organize variables: Australia
#######################

## wave 1

# Rename redistribution items to intuitive names

au1$redist1 <- au1$Q4_1
au1$redist2 <- au1$Q4_2
au1$redist3 <- au1$Q4_3
au1$redist4 <- au1$Q4_4
au1$redist5 <- au1$Q4_5
au1$redist6 <- au1$Q4_6

# Recode redistribution variables to numeric.

au1$redist1 <- as.numeric(mapvalues(au1$redist1,c("Strongly disagree","Somewhat disagree","Neither agree nor disagree",
                                                  "Somewhat agree","Strongly agree"),
                                    1:5))
au1$redist2 <- as.numeric(mapvalues(au1$redist2,c("Strongly disagree","Somewhat disagree","Neither agree nor disagree",
                                                  "Somewhat agree","Strongly agree"),
                                    5:1))
au1$redist3 <- as.numeric(mapvalues(au1$redist3,c("Strongly disagree","Somewhat disagree","Neither agree nor disagree",
                                                  "Somewhat agree","Strongly agree"),
                                    5:1))
au1$redist4 <- as.numeric(mapvalues(au1$redist4,c("Strongly disagree","Somewhat disagree","Neither agree nor disagree",
                                                  "Somewhat agree","Strongly agree"),
                                    5:1))
au1$redist5 <- as.numeric(mapvalues(au1$redist5,c("Strongly disagree","Somewhat disagree","Neither agree nor disagree",
                                                  "Somewhat agree","Strongly agree"),
                                    1:5))
au1$redist6 <- as.numeric(mapvalues(au1$redist6,c("Strongly disagree","Somewhat disagree","Neither agree nor disagree",
                                                  "Somewhat agree","Strongly agree"),
                                    1:5))

# Create average measure of redistribution attitudes

au1$redist_avg_w1 <- apply(au1[,paste("redist",1:6,sep="")],1,mean,na.rm=T)



# Form binary variable of news consumption

au1$news_frequency <- ifelse(au1$Q9_1s == "6 days" | 
                               au1$Q9_1s == "7 days" |
                               au1$Q9_2s == "6 days" | 
                               au1$Q9_2s == "7 days" |
                               au1$Q9_3s == "6 days" | 
                               au1$Q9_3s == "7 days" |
                               au1$Q9_4s == "6 days" | 
                               au1$Q9_4s == "7 days" |
                               au1$Q9_5s == "6 days" | 
                               au1$Q9_5s == "7 days", 1, 0)

au1$news_frequency[is.na(au1$news_frequency)] <- 0
au1$news_frequency = factor(au1$news_frequency, levels = c(0, 1), 
                            labels = c('Low Exposure to Financial News', 'High Exposure to Financial News'))



# Produce mean and factor score of racial resentment variable

au1$rr_avg <- apply(au1[,paste("rr",1:4,sep="")],1,mean,na.rm=T)

au1$rr_avg_polychoric_factor <- fa(au1[,paste(c("rr1","rr2", "rr3", "rr4"))],nfactors=1,rotate="Promax", cor = "poly")$scores[,1]
au1$rr_avg_polychoric_factor <- scales:::rescale(au1$rr_avg_polychoric_factor,to=c(1,5))


# Produce mean and factor score of individualism scale items

au1$indiv_avg <- apply(au1[,paste("indiv",1:5,sep="")],1,mean,na.rm=T)

au1$indiv_avg_polychoric_factor <- fa(au1[,paste(c("indiv1","indiv2", "indiv3", "indiv4", "indiv5"))],nfactors=1,rotate="Promax", cor = "poly")$scores[,1]
au1$indiv_avg_polychoric_factor <- scales:::rescale(au1$indiv_avg_polychoric_factor,to=c(1,5))



# Produce affective polarization / Partisan identification construct

au1$thermo_lab <- au1$Q30a_g_1 #Labour
au1$thermo_lib <- au1$Q30a_g_2 #Liberal
au1$thermo_lab <- as.numeric(mapvalues(au1$thermo_lab, c("Extremely cold0","1","2","3","4",
                                                         "Neither warm nor cold5","6","7","8",
                                                         "9","Extremely warm10"), 0:10))
au1$thermo_lib <- as.numeric(mapvalues(au1$thermo_lib, c("Extremely cold0","1","2","3","4",
                                                         "Neither warm nor cold5","6","7","8",
                                                         "9","Extremely warm10"), 0:10))
au1$thermo_diff <- au1$thermo_lib - au1$thermo_lab


## wave 3

# Clean redistribution variable - change item names, recode to numeric, produce mean

au3$redist1 <- au3$Q4_1
au3$redist2 <- au3$Q4_2
au3$redist3 <- au3$Q4_3
au3$redist4 <- au3$Q4_4
au3$redist5 <- au3$Q4_5
au3$redist6 <- au3$Q4_6
au3$redist1 <- as.numeric(mapvalues(au3$redist1,c("Strongly disagree","Somewhat disagree","Neither agree nor disagree",
                                                  "Somewhat agree","Strongly agree"),
                                    1:5))
au3$redist2 <- as.numeric(mapvalues(au3$redist2,c("Strongly disagree","Somewhat disagree","Neither agree nor disagree",
                                                  "Somewhat agree","Strongly agree"),
                                    1:5))
au3$redist3 <- as.numeric(mapvalues(au3$redist3,c("Strongly disagree","Somewhat disagree","Neither agree nor disagree",
                                                  "Somewhat agree","Strongly agree"),
                                    1:5))
au3$redist4 <- as.numeric(mapvalues(au3$redist4,c("Strongly disagree","Somewhat disagree","Neither agree nor disagree",
                                                  "Somewhat agree","Strongly agree"),
                                    5:1))
au3$redist5 <- as.numeric(mapvalues(au3$redist5,c("Strongly disagree","Somewhat disagree","Neither agree nor disagree",
                                                  "Somewhat agree","Strongly agree"),
                                    1:5))
au3$redist6 <- as.numeric(mapvalues(au3$redist6,c("Strongly disagree","Somewhat disagree","Neither agree nor disagree",
                                                  "Somewhat agree","Strongly agree"),
                                    1:5))
au3$redist_avg <- apply(au3[,paste("redist",1:6,sep="")],1,mean,na.rm=T)

au3$redist_fac <- fa(au3[,paste("redist",1:6,sep="")],nfactors=1,rotate="Promax", cor = "poly")$scores[,1]
au3$redist_fac_1to5 <- scales:::rescale(au3$redist_fac,to=c(1,5))

# Clean elite domination variable - change item names, recode to numeric, produce mean and factor score

au3$rigged1 <- au3$Q19_1
au3$rigged2 <- au3$Q19_2
au3$rigged3 <- au3$Q19_3
au3$rigged4 <- au3$Q19_4
au3$rigged5 <- au3$Q19_5
au3$rigged6 <- au3$Q19_6
au3$rigged1 <- as.numeric(mapvalues(au3$rigged1,c("Strongly agree","Somewhat agree",
                                                  "Neither agree nor disagree","Somewhat disagree",
                                                  "Strongly disagree"), 5:1))
au3$rigged2 <- as.numeric(mapvalues(au3$rigged2,c("Strongly agree","Somewhat agree",
                                                  "Neither agree nor disagree","Somewhat disagree",
                                                  "Strongly disagree"), 5:1))
au3$rigged3 <- as.numeric(mapvalues(au3$rigged3,c("Strongly agree","Somewhat agree",
                                                  "Neither agree nor disagree","Somewhat disagree",
                                                  "Strongly disagree"), 5:1))
au3$rigged4 <- as.numeric(mapvalues(au3$rigged4,c("Strongly agree","Somewhat agree",
                                                  "Neither agree nor disagree","Somewhat disagree",
                                                  "Strongly disagree"), 5:1))
au3$rigged5 <- as.numeric(mapvalues(au3$rigged5,c("Strongly agree","Somewhat agree",
                                                  "Neither agree nor disagree","Somewhat disagree",
                                                  "Strongly disagree"), 1:5))
au3$rigged6 <- as.numeric(mapvalues(au3$rigged6,c("Strongly agree","Somewhat agree",
                                                  "Neither agree nor disagree","Somewhat disagree",
                                                  "Strongly disagree"), 5:1))

au3$elitedom <- apply(au3[,paste("rigged",1:6,sep="")],1,mean,na.rm=T)

au3$elitedom_polychoric_factor_score <- fa(au3[,paste("rigged",1:6,sep="")],nfactors=1,rotate="Promax", cor = "poly")$scores[,1]
au3$elitedom_polychoric_factor_score <- scales:::rescale(au3$elitedom_polychoric_factor_score,to=c(1,5))









#######################
#### Organize variables: France
#######################

## wave 1

# Rename redistribution items to intuitive names

fr1$redist1 <- fr1$Q4_1
fr1$redist2 <- fr1$Q4_2
fr1$redist3 <- fr1$Q4_3
fr1$redist4 <- fr1$Q4_4
fr1$redist5 <- fr1$Q4_5
fr1$redist6 <- fr1$Q4_6

# Recode redistribution variables to numeric.

fr1$redist1 <- as.numeric(mapvalues(fr1$redist1,c("Pas du tout d'accord","Pas vraiment d'accord","Ni d'accord, ni pas d'accord",
                                                  "Plutôt d'accord","Tout à fait d'accord"),
                                    1:5))
fr1$redist2 <- as.numeric(mapvalues(fr1$redist2,c("Pas du tout d'accord","Pas vraiment d'accord","Ni d'accord, ni pas d'accord",
                                                  "Plutôt d'accord","Tout à fait d'accord"),
                                    5:1))
fr1$redist3 <- as.numeric(mapvalues(fr1$redist3,c("Pas du tout d'accord","Pas vraiment d'accord","Ni d'accord, ni pas d'accord",
                                                  "Plutôt d'accord","Tout à fait d'accord"),
                                    5:1))
fr1$redist4 <- as.numeric(mapvalues(fr1$redist4,c("Pas du tout d'accord","Pas vraiment d'accord","Ni d'accord, ni pas d'accord",
                                                  "Plutôt d'accord","Tout à fait d'accord"),
                                    5:1))
fr1$redist5 <- as.numeric(mapvalues(fr1$redist5,c("Pas du tout d'accord","Pas vraiment d'accord","Ni d'accord, ni pas d'accord",
                                                  "Plutôt d'accord","Tout à fait d'accord"),
                                    1:5))
fr1$redist6 <- as.numeric(mapvalues(fr1$redist6,c("Pas du tout d'accord","Pas vraiment d'accord","Ni d'accord, ni pas d'accord",
                                                  "Plutôt d'accord","Tout à fait d'accord"),
                                    1:5))

# Create average measure of redistribution attitudes

fr1$redist_avg_w1 <- apply(fr1[,paste("redist",1:6,sep="")],1,mean,na.rm=T)



# Form binary variable of news consumption

fr1$news_frequency <- ifelse(fr1$Q9_1s == "6 jours" | 
                               fr1$Q9_1s == "7 jours" |
                               fr1$Q9_2s == "6 jours" | 
                               fr1$Q9_2s == "7 jours" |
                               fr1$Q9_3s == "6 jours" | 
                               fr1$Q9_3s == "7 jours" |
                               fr1$Q9_4s == "6 jours" | 
                               fr1$Q9_4s == "7 jours" |
                               fr1$Q9_5s == "6 jours" | 
                               fr1$Q9_5s == "7 jours", 1, 0)

fr1$news_frequency[is.na(fr1$news_frequency)] <- 0
fr1$news_frequency = factor(fr1$news_frequency, levels = c(0, 1), 
                            labels = c('Low Exposure to Financial News', 'High Exposure to Financial News'))



# Produce mean and factor score of racial resentment variable

fr1$rr_avg <- apply(fr1[,paste("rr",1:4,sep="")],1,mean,na.rm=T)

fr1$rr_avg_polychoric_factor <- fa(fr1[,paste(c("rr1","rr2", "rr3", "rr4"))],nfactors=1,rotate="Promax", cor = "poly")$scores[,1]
fr1$rr_avg_polychoric_factor <- scales:::rescale(fr1$rr_avg_polychoric_factor,to=c(1,5))


# Produce mean and factor score of individualism scale items

fr1$indiv_avg <- apply(fr1[,paste("indiv",1:5,sep="")],1,mean,na.rm=T)

fr1$indiv_avg_polychoric_factor <- fa(fr1[,paste(c("indiv1","indiv2", "indiv3", "indiv4", "indiv5"))],nfactors=1,rotate="Promax", cor = "poly")$scores[,1]
fr1$indiv_avg_polychoric_factor <- scales:::rescale(fr1$indiv_avg_polychoric_factor,to=c(1,5))


# Produce affective polarization / partisan identification construct 

fr1$thermo_lab <- fr1$Q30a_g_1 #En Marche
fr1$thermo_lib <- fr1$Q30a_g_2 #Les Republicains
fr1$thermo_lab <- as.numeric(mapvalues(fr1$thermo_lab, c("Extrêmement négative0","1","2","3","4",
                                                         "Ni positive, ni négative5","6","7","8",
                                                         "9","Extrêmement positive10"), 0:10))
fr1$thermo_lib <- as.numeric(mapvalues(fr1$thermo_lib, c("Extrêmement négative0","1","2","3","4",
                                                         "Ni positive, ni négative5","6","7","8",
                                                         "9","Extrêmement positive10"), 0:10))
fr1$thermo_diff <- fr1$thermo_lib - fr1$thermo_lab



## wave 3

# Clean redistribution variable - change item names, recode to numeric, produce mean

fr3$redist1 <- fr3$Q4_1
fr3$redist2 <- fr3$Q4_2
fr3$redist3 <- fr3$Q4_3
fr3$redist4 <- fr3$Q4_4
fr3$redist5 <- fr3$Q4_5
fr3$redist6 <- fr3$Q4_6
fr3$redist1 <- as.numeric(mapvalues(fr3$redist1,c("Pas du tout d'accord","Pas vraiment d'accord","Ni d'accord, ni pas d'accord",
                                                  "Plutôt d'accord","Tout à fait d'accord"),
                                    1:5))
fr3$redist2 <- as.numeric(mapvalues(fr3$redist2,c("Pas du tout d'accord","Pas vraiment d'accord","Ni d'accord, ni pas d'accord",
                                                  "Plutôt d'accord","Tout à fait d'accord"),
                                    1:5))
fr3$redist3 <- as.numeric(mapvalues(fr3$redist3,c("Pas du tout d'accord","Pas vraiment d'accord","Ni d'accord, ni pas d'accord",
                                                  "Plutôt d'accord","Tout à fait d'accord"),
                                    1:5))
fr3$redist4 <- as.numeric(mapvalues(fr3$redist4,c("Pas du tout d'accord","Pas vraiment d'accord","Ni d'accord, ni pas d'accord",
                                                  "Plutôt d'accord","Tout à fait d'accord"),
                                    5:1))
fr3$redist5 <- as.numeric(mapvalues(fr3$redist5,c("Pas du tout d'accord","Pas vraiment d'accord","Ni d'accord, ni pas d'accord",
                                                  "Plutôt d'accord","Tout à fait d'accord"),
                                    1:5))
fr3$redist6 <- as.numeric(mapvalues(fr3$redist6,c("Pas du tout d'accord","Pas vraiment d'accord","Ni d'accord, ni pas d'accord",
                                                  "Plutôt d'accord","Tout à fait d'accord"),
                                    1:5))
fr3$redist_avg <- apply(fr3[,paste("redist",1:6,sep="")],1,mean,na.rm=T)

fr3$redist_fac <- fa(fr3[,paste("redist",1:6,sep="")],nfactors=1,rotate="Promax", cor = "poly")$scores[,1]
fr3$redist_fac_1to5 <- scales:::rescale(fr3$redist_fac,to=c(1,5))


# Clean elite domination variable - change item names, recode to numeric, produce mean and factor score

fr3$rigged1 <- fr3$Q19_1
fr3$rigged2 <- fr3$Q19_2
fr3$rigged3 <- fr3$Q19_3
fr3$rigged4 <- fr3$Q19_4
fr3$rigged5 <- fr3$Q19_5
fr3$rigged6 <- fr3$Q19_6
fr3$rigged1 <- as.numeric(mapvalues(fr3$rigged1,c("Tout à fait d'accord","Plutôt d'accord",
                                                  "Ni d'accord, ni pas d'accord","Pas vraiment d'accord",
                                                  "Pas du tout d'accord"), 5:1))
fr3$rigged2 <- as.numeric(mapvalues(fr3$rigged2,c("Tout à fait d'accord","Plutôt d'accord",
                                                  "Ni d'accord, ni pas d'accord","Pas vraiment d'accord",
                                                  "Pas du tout d'accord"), 5:1))
fr3$rigged3 <- as.numeric(mapvalues(fr3$rigged3,c("Tout à fait d'accord","Plutôt d'accord",
                                                  "Ni d'accord, ni pas d'accord","Pas vraiment d'accord",
                                                  "Pas du tout d'accord"), 5:1))
fr3$rigged4 <- as.numeric(mapvalues(fr3$rigged4,c("Tout à fait d'accord","Plutôt d'accord",
                                                  "Ni d'accord, ni pas d'accord","Pas vraiment d'accord",
                                                  "Pas du tout d'accord"), 5:1))
fr3$rigged5 <- as.numeric(mapvalues(fr3$rigged5,c("Tout à fait d'accord","Plutôt d'accord",
                                                  "Ni d'accord, ni pas d'accord","Pas vraiment d'accord",
                                                  "Pas du tout d'accord"), 1:5))
fr3$rigged6 <- as.numeric(mapvalues(fr3$rigged6,c("Tout à fait d'accord","Plutôt d'accord",
                                                  "Ni d'accord, ni pas d'accord","Pas vraiment d'accord",
                                                  "Pas du tout d'accord"), 5:1))

fr3$elitedom <- apply(fr3[,paste("rigged",1:6,sep="")],1,mean,na.rm=T)

fr3$elitedom_polychoric_factor_score <- fa(fr3[,paste("rigged",1:6,sep="")],nfactors=1,rotate="Promax", cor = "poly")$scores[,1]
fr3$elitedom_polychoric_factor_score <- scales:::rescale(fr3$elitedom_polychoric_factor_score,to=c(1,5))









#######################
#### Organize variables: Germany
#######################

## wave 1

# Rename redistribution items to intuitive names

de1$redist1 <- de1$Q4_1
de1$redist2 <- de1$Q4_2
de1$redist3 <- de1$Q4_3
de1$redist4 <- de1$Q4_4
de1$redist5 <- de1$Q4_5
de1$redist6 <- de1$Q4_6

# Recode redistribution variables to numeric.

de1$redist1 <- as.numeric(mapvalues(de1$redist1,c("Stimme überhaupt nicht zu","Stimme eher nicht zu","Neutral",
                                                  "Stimme eher zu","Stimme vollkommen zu"),
                                    1:5))
de1$redist2 <- as.numeric(mapvalues(de1$redist2,c("Stimme überhaupt nicht zu","Stimme eher nicht zu","Neutral",
                                                  "Stimme eher zu","Stimme vollkommen zu"),
                                    5:1))
de1$redist3 <- as.numeric(mapvalues(de1$redist3,c("Stimme überhaupt nicht zu","Stimme eher nicht zu","Neutral",
                                                  "Stimme eher zu","Stimme vollkommen zu"),
                                    5:1))
de1$redist4 <- as.numeric(mapvalues(de1$redist4,c("Stimme überhaupt nicht zu","Stimme eher nicht zu","Neutral",
                                                  "Stimme eher zu","Stimme vollkommen zu"),
                                    5:1))
de1$redist5 <- as.numeric(mapvalues(de1$redist5,c("Stimme überhaupt nicht zu","Stimme eher nicht zu","Neutral",
                                                  "Stimme eher zu","Stimme vollkommen zu"),
                                    1:5))
de1$redist6 <- as.numeric(mapvalues(de1$redist6,c("Stimme überhaupt nicht zu","Stimme eher nicht zu","Neutral",
                                                  "Stimme eher zu","Stimme vollkommen zu"),
                                    1:5))

# Create average measure of redistribution attitudes

de1$redist_avg_w1 <- apply(de1[,paste("redist",1:6,sep="")],1,mean,na.rm=T)



# Form binary variable of news consumption

de1$news_frequency <- ifelse(de1$Q9_1s == "An 6 Tagen" | 
                               de1$Q9_1s == "An 7 Tagen" |
                               de1$Q9_2s == "An 6 Tagen" | 
                               de1$Q9_2s == "An 7 Tagen" |
                               de1$Q9_3s == "An 6 Tagen" | 
                               de1$Q9_3s == "An 7 Tagen" |
                               de1$Q9_4s == "An 6 Tagen" | 
                               de1$Q9_4s == "An 7 Tagen" |
                               de1$Q9_5s == "An 6 Tagen" | 
                               de1$Q9_5s == "An 7 Tagen", 1, 0)

de1$news_frequency[is.na(de1$news_frequency)] <- 0
de1$news_frequency = factor(de1$news_frequency, levels = c(0, 1), 
                            labels = c('Low Exposure to Financial News', 'High Exposure to Financial News'))



# Produce mean and factor score of racial resentment variable

de1$rr_avg <- apply(de1[,paste("rr",1:4,sep="")],1,mean,na.rm=T)

de1$rr_avg_polychoric_factor <- fa(de1[,paste(c("rr1","rr2", "rr3", "rr4"))],nfactors=1,rotate="Promax", cor = "poly")$scores[,1]
de1$rr_avg_polychoric_factor <- scales:::rescale(de1$rr_avg_polychoric_factor,to=c(1,5))


# Produce mean and factor score of individualism scale items

de1$indiv_avg <- apply(de1[,paste("indiv",1:5,sep="")],1,mean,na.rm=T)

de1$indiv_avg_polychoric_factor <- fa(de1[,paste(c("indiv1","indiv2", "indiv3", "indiv4", "indiv5"))],nfactors=1,rotate="Promax", cor = "poly")$scores[,1]
de1$indiv_avg_polychoric_factor <- scales:::rescale(de1$indiv_avg_polychoric_factor,to=c(1,5))


# Produce affective polarization / partisan identification construct 

de1$thermo_lab <- de1$Q30a_g_2 #SPD
de1$thermo_lib <- de1$Q30a_g_1 #CDU
de1$thermo_lab <- as.numeric(mapvalues(de1$thermo_lab, c("Sehr stark abgeneigt0","1","2","3","4",
                                                         "Neutral5","6","7","8",
                                                         "9","Sehr stark zugeneigt10"), 0:10))
de1$thermo_lib <- as.numeric(mapvalues(de1$thermo_lib, c("Sehr stark abgeneigt0","1","2","3","4",
                                                         "Neutral5","6","7","8",
                                                         "9","Sehr stark zugeneigt10"), 0:10))
de1$thermo_diff <- de1$thermo_lib - de1$thermo_lab


## wave 3

# Clean redistribution variable - change item names, recode to numeric, produce mean

de3$redist1 <- de3$Q4_1
de3$redist2 <- de3$Q4_2
de3$redist3 <- de3$Q4_3
de3$redist4 <- de3$Q4_4
de3$redist5 <- de3$Q4_5
de3$redist6 <- de3$Q4_6
de3$redist1 <- as.numeric(mapvalues(de3$redist1,c("Stimme überhaupt nicht zu","Stimme eher nicht zu","Neutral",
                                                  "Stimme eher zu","Stimme vollkommen zu"),
                                    1:5))
de3$redist2 <- as.numeric(mapvalues(de3$redist2,c("Stimme überhaupt nicht zu","Stimme eher nicht zu","Neutral",
                                                  "Stimme eher zu","Stimme vollkommen zu"),
                                    1:5))
de3$redist3 <- as.numeric(mapvalues(de3$redist3,c("Stimme überhaupt nicht zu","Stimme eher nicht zu","Neutral",
                                                  "Stimme eher zu","Stimme vollkommen zu"),
                                    1:5))
de3$redist4 <- as.numeric(mapvalues(de3$redist4,c("Stimme überhaupt nicht zu","Stimme eher nicht zu","Neutral",
                                                  "Stimme eher zu","Stimme vollkommen zu"),
                                    5:1))
de3$redist5 <- as.numeric(mapvalues(de3$redist5,c("Stimme überhaupt nicht zu","Stimme eher nicht zu","Neutral",
                                                  "Stimme eher zu","Stimme vollkommen zu"),
                                    1:5))
de3$redist6 <- as.numeric(mapvalues(de3$redist6,c("Stimme überhaupt nicht zu","Stimme eher nicht zu","Neutral",
                                                  "Stimme eher zu","Stimme vollkommen zu"),
                                    1:5))
de3$redist_avg <- apply(de3[,paste("redist",1:6,sep="")],1,mean,na.rm=T)

de3$redist_fac <- fa(de3[,paste("redist",1:6,sep="")],nfactors=1,rotate="Promax", cor = "poly")$scores[,1]
de3$redist_fac_1to5 <- scales:::rescale(de3$redist_fac,to=c(1,5))


# Clean elite domination variable - change item names, recode to numeric, produce mean and factor score

de3$rigged1 <- de3$Q19_1
de3$rigged2 <- de3$Q19_2
de3$rigged3 <- de3$Q19_3
de3$rigged4 <- de3$Q19_4
de3$rigged5 <- de3$Q19_5
de3$rigged6 <- de3$Q19_6
de3$rigged1 <- as.numeric(mapvalues(de3$rigged1,c("Stimme vollkommen zu","Stimme eher zu",
                                                  "Neutral","Stimme eher nicht zu",
                                                  "Stimme überhaupt nicht zu"), 5:1))
de3$rigged2 <- as.numeric(mapvalues(de3$rigged2,c("Stimme vollkommen zu","Stimme eher zu",
                                                  "Neutral","Stimme eher nicht zu",
                                                  "Stimme überhaupt nicht zu"), 5:1))
de3$rigged3 <- as.numeric(mapvalues(de3$rigged3,c("Stimme vollkommen zu","Stimme eher zu",
                                                  "Neutral","Stimme eher nicht zu",
                                                  "Stimme überhaupt nicht zu"), 5:1))
de3$rigged4 <- as.numeric(mapvalues(de3$rigged4,c("Stimme vollkommen zu","Stimme eher zu",
                                                  "Neutral","Stimme eher nicht zu",
                                                  "Stimme überhaupt nicht zu"), 5:1))
de3$rigged5 <- as.numeric(mapvalues(de3$rigged5,c("Stimme vollkommen zu","Stimme eher zu",
                                                  "Neutral","Stimme eher nicht zu",
                                                  "Stimme überhaupt nicht zu"), 1:5))
de3$rigged6 <- as.numeric(mapvalues(de3$rigged6,c("Stimme vollkommen zu","Stimme eher zu",
                                                  "Neutral","Stimme eher nicht zu",
                                                  "Stimme überhaupt nicht zu"), 5:1))

de3$elitedom <- apply(de3[,paste("rigged",1:6,sep="")],1,mean,na.rm=T)

de3$elitedom_polychoric_factor_score <- fa(de3[,paste("rigged",1:6,sep="")],nfactors=1,rotate="Promax", cor = "poly")$scores[,1]
de3$elitedom_polychoric_factor_score <- scales:::rescale(de3$elitedom_polychoric_factor_score,to=c(1,5))






#######################
#### Organize variables: United Kingdom
#######################

## wave 1

# Rename redistribution items to intuitive names

uk1$redist1 <- uk1$Q4_1
uk1$redist2 <- uk1$Q4_2
uk1$redist3 <- uk1$Q4_3
uk1$redist4 <- uk1$Q4_4
uk1$redist5 <- uk1$Q4_5
uk1$redist6 <- uk1$Q4_6

# Recode redistribution variables to numeric.

uk1$redist1 <- as.numeric(mapvalues(uk1$redist1,c("Strongly disagree","Somewhat disagree","Neither agree nor disagree",
                                                  "Somewhat agree","Strongly agree"),
                                    1:5))
uk1$redist2 <- as.numeric(mapvalues(uk1$redist2,c("Strongly disagree","Somewhat disagree","Neither agree nor disagree",
                                                  "Somewhat agree","Strongly agree"),
                                    5:1))
uk1$redist3 <- as.numeric(mapvalues(uk1$redist3,c("Strongly disagree","Somewhat disagree","Neither agree nor disagree",
                                                  "Somewhat agree","Strongly agree"),
                                    5:1))
uk1$redist4 <- as.numeric(mapvalues(uk1$redist4,c("Strongly disagree","Somewhat disagree","Neither agree nor disagree",
                                                  "Somewhat agree","Strongly agree"),
                                    5:1))
uk1$redist5 <- as.numeric(mapvalues(uk1$redist5,c("Strongly disagree","Somewhat disagree","Neither agree nor disagree",
                                                  "Somewhat agree","Strongly agree"),
                                    1:5))
uk1$redist6 <- as.numeric(mapvalues(uk1$redist6,c("Strongly disagree","Somewhat disagree","Neither agree nor disagree",
                                                  "Somewhat agree","Strongly agree"),
                                    1:5))

# Create average measure of redistribution attitudes

uk1$redist_avg_w1 <- apply(uk1[,paste("redist",1:6,sep="")],1,mean,na.rm=T)



# Form binary variable of news consumption

uk1$news_frequency <- ifelse(uk1$Q9_1s == "6 days" | 
                               uk1$Q9_1s == "7 days" |
                               uk1$Q9_2s == "6 days" | 
                               uk1$Q9_2s == "7 days" |
                               uk1$Q9_3s == "6 days" | 
                               uk1$Q9_3s == "7 days" |
                               uk1$Q9_4s == "6 days" | 
                               uk1$Q9_4s == "7 days" |
                               uk1$Q9_5s == "6 days" | 
                               uk1$Q9_5s == "7 days", 1, 0)

uk1$news_frequency[is.na(uk1$news_frequency)] <- 0
uk1$news_frequency = factor(uk1$news_frequency, levels = c(0, 1), 
                            labels = c('Low Exposure to Financial News', 'High Exposure to Financial News'))



# Produce mean and factor score of racial resentment variable

uk1$rr_avg <- apply(uk1[,paste("rr",1:4,sep="")],1,mean,na.rm=T)

uk1$rr_avg_polychoric_factor <- fa(uk1[,paste(c("rr1","rr2", "rr3", "rr4"))],nfactors=1,rotate="Promax", cor = "poly")$scores[,1]
uk1$rr_avg_polychoric_factor <- scales:::rescale(uk1$rr_avg_polychoric_factor,to=c(1,5))


# Produce mean and factor score of individualism scale items

uk1$indiv_avg <- apply(uk1[,paste("indiv",1:5,sep="")],1,mean,na.rm=T)

uk1$indiv_avg_polychoric_factor <- fa(uk1[,paste(c("indiv1","indiv2", "indiv3", "indiv4", "indiv5"))],nfactors=1,rotate="Promax", cor = "poly")$scores[,1]
uk1$indiv_avg_polychoric_factor <- scales:::rescale(uk1$indiv_avg_polychoric_factor,to=c(1,5))


# Produce affective polarization / partisan identification construct 

uk1$thermo_lab <- uk1$Q30a_g_1 #Labour
uk1$thermo_lib <- uk1$Q30a_g_2 #Conservatives
uk1$thermo_lab <- as.numeric(mapvalues(uk1$thermo_lab, c("Extremely cold0","1","2","3","4",
                                                         "Neither warm nor cold5","6","7","8",
                                                         "9","Extremely warm10"), 0:10))
uk1$thermo_lib <- as.numeric(mapvalues(uk1$thermo_lib, c("Extremely cold0","1","2","3","4",
                                                         "Neither warm nor cold5","6","7","8",
                                                         "9","Extremely warm10"), 0:10))
uk1$thermo_diff <- uk1$thermo_lib - uk1$thermo_lab


## wave 3

# Clean redistribution variable - change item names, recode to numeric, produce mean

uk3$redist1 <- uk3$Q4_1
uk3$redist2 <- uk3$Q4_2
uk3$redist3 <- uk3$Q4_3
uk3$redist4 <- uk3$Q4_4
uk3$redist5 <- uk3$Q4_5
uk3$redist6 <- uk3$Q4_6
uk3$redist1 <- as.numeric(mapvalues(uk3$redist1,c("Strongly disagree","Somewhat disagree","Neither agree nor disagree",
                                                  "Somewhat agree","Strongly agree"),
                                    1:5))
uk3$redist2 <- as.numeric(mapvalues(uk3$redist2,c("Strongly disagree","Somewhat disagree","Neither agree nor disagree",
                                                  "Somewhat agree","Strongly agree"),
                                    1:5))
uk3$redist3 <- as.numeric(mapvalues(uk3$redist3,c("Strongly disagree","Somewhat disagree","Neither agree nor disagree",
                                                  "Somewhat agree","Strongly agree"),
                                    1:5))
uk3$redist4 <- as.numeric(mapvalues(uk3$redist4,c("Strongly disagree","Somewhat disagree","Neither agree nor disagree",
                                                  "Somewhat agree","Strongly agree"),
                                    5:1))
uk3$redist5 <- as.numeric(mapvalues(uk3$redist5,c("Strongly disagree","Somewhat disagree","Neither agree nor disagree",
                                                  "Somewhat agree","Strongly agree"),
                                    1:5))
uk3$redist6 <- as.numeric(mapvalues(uk3$redist6,c("Strongly disagree","Somewhat disagree","Neither agree nor disagree",
                                                  "Somewhat agree","Strongly agree"),
                                    1:5))
uk3$redist_avg <- apply(uk3[,paste("redist",1:6,sep="")],1,mean,na.rm=T)

uk3$redist_fac <- fa(uk3[,paste("redist",1:6,sep="")],nfactors=1,rotate="Promax", cor = "poly")$scores[,1]
uk3$redist_fac_1to5 <- scales:::rescale(uk3$redist_fac,to=c(1,5))


# Clean elite domination variable - change item names, recode to numeric, produce mean and factor score

uk3$rigged1 <- uk3$Q19_1
uk3$rigged2 <- uk3$Q19_2
uk3$rigged3 <- uk3$Q19_3
uk3$rigged4 <- uk3$Q19_4
uk3$rigged5 <- uk3$Q19_5
uk3$rigged6 <- uk3$Q19_6
uk3$rigged1 <- as.numeric(mapvalues(uk3$rigged1,c("Strongly agree","Somewhat agree",
                                                  "Neither agree nor disagree","Somewhat disagree",
                                                  "Strongly disagree"), 5:1))
uk3$rigged2 <- as.numeric(mapvalues(uk3$rigged2,c("Strongly agree","Somewhat agree",
                                                  "Neither agree nor disagree","Somewhat disagree",
                                                  "Strongly disagree"), 5:1))
uk3$rigged3 <- as.numeric(mapvalues(uk3$rigged3,c("Strongly agree","Somewhat agree",
                                                  "Neither agree nor disagree","Somewhat disagree",
                                                  "Strongly disagree"), 5:1))
uk3$rigged4 <- as.numeric(mapvalues(uk3$rigged4,c("Strongly agree","Somewhat agree",
                                                  "Neither agree nor disagree","Somewhat disagree",
                                                  "Strongly disagree"), 5:1))
uk3$rigged5 <- as.numeric(mapvalues(uk3$rigged5,c("Strongly agree","Somewhat agree",
                                                  "Neither agree nor disagree","Somewhat disagree",
                                                  "Strongly disagree"), 1:5))
uk3$rigged6 <- as.numeric(mapvalues(uk3$rigged6,c("Strongly agree","Somewhat agree",
                                                  "Neither agree nor disagree","Somewhat disagree",
                                                  "Strongly disagree"), 5:1))

uk3$elitedom <- apply(uk3[,paste("rigged",1:6,sep="")],1,mean,na.rm=T)

uk3$elitedom_polychoric_factor_score <- fa(uk3[,paste("rigged",1:6,sep="")],nfactors=1,rotate="Promax", cor = "poly")$scores[,1]
uk3$elitedom_polychoric_factor_score <- scales:::rescale(uk3$elitedom_polychoric_factor_score,to=c(1,5))







#######################
#### Organize variables: United States
#######################

## wave 1

# Rename redistribution items to intuitive names

us1$redist1 <- us1$Q4_1
us1$redist2 <- us1$Q4_2
us1$redist3 <- us1$Q4_3
us1$redist4 <- us1$Q4_4
us1$redist5 <- us1$Q4_5
us1$redist6 <- us1$Q4_6

# Recode redistribution variables to numeric.

us1$redist1 <- as.numeric(mapvalues(us1$redist1,c("Strongly disagree","Somewhat disagree","Neither agree nor disagree",
                                                  "Somewhat agree","Strongly agree"),
                                    1:5))
us1$redist2 <- as.numeric(mapvalues(us1$redist2,c("Strongly disagree","Somewhat disagree","Neither agree nor disagree",
                                                  "Somewhat agree","Strongly agree"),
                                    5:1))
us1$redist3 <- as.numeric(mapvalues(us1$redist3,c("Strongly disagree","Somewhat disagree","Neither agree nor disagree",
                                                  "Somewhat agree","Strongly agree"),
                                    5:1))
us1$redist4 <- as.numeric(mapvalues(us1$redist4,c("Strongly disagree","Somewhat disagree","Neither agree nor disagree",
                                                  "Somewhat agree","Strongly agree"),
                                    5:1))
us1$redist5 <- as.numeric(mapvalues(us1$redist5,c("Strongly disagree","Somewhat disagree","Neither agree nor disagree",
                                                  "Somewhat agree","Strongly agree"),
                                    1:5))
us1$redist6 <- as.numeric(mapvalues(us1$redist6,c("Strongly disagree","Somewhat disagree","Neither agree nor disagree",
                                                  "Somewhat agree","Strongly agree"),
                                    1:5))

# Create average measure of redistribution attitudes

us1$redist_avg_w1 <- apply(us1[,paste("redist",1:6,sep="")],1,mean,na.rm=T)



# Form binary variable of news consumption

us1$news_frequency <- ifelse(us1$Q9_1s == "6 days" | 
                               us1$Q9_1s == "7 days" |
                               us1$Q9_2s == "6 days" | 
                               us1$Q9_2s == "7 days" |
                               us1$Q9_3s == "6 days" | 
                               us1$Q9_3s == "7 days" |
                               us1$Q9_4s == "6 days" | 
                               us1$Q9_4s == "7 days" |
                               us1$Q9_5s == "6 days" | 
                               us1$Q9_5s == "7 days", 1, 0)

us1$news_frequency[is.na(us1$news_frequency)] <- 0
us1$news_frequency = factor(us1$news_frequency, levels = c(0, 1), 
                            labels = c('Low Exposure to Financial News', 'High Exposure to Financial News'))



# Produce mean and factor score of racial resentment variable

us1$rr_avg <- apply(us1[,paste("rr",1:4,sep="")],1,mean,na.rm=T)

us1$rr_avg_polychoric_factor <- fa(us1[,paste(c("rr1","rr2", "rr3", "rr4"))],nfactors=1,rotate="Promax", cor = "poly")$scores[,1]
us1$rr_avg_polychoric_factor <- scales:::rescale(us1$rr_avg_polychoric_factor,to=c(1,5))


# Produce mean and factor score of individualism scale items

us1$indiv_avg <- apply(us1[,paste("indiv",1:5,sep="")],1,mean,na.rm=T)

us1$indiv_avg_polychoric_factor <- fa(us1[,paste(c("indiv1","indiv2", "indiv3", "indiv4", "indiv5"))],nfactors=1,rotate="Promax", cor = "poly")$scores[,1]
us1$indiv_avg_polychoric_factor <- scales:::rescale(us1$indiv_avg_polychoric_factor,to=c(1,5))


# Produce affective polarization / partisan identification construct 

us1$thermo_lab <- us1$Q30a_g_1 #Dems
us1$thermo_lib <- us1$Q30a_g_2 #GOP
us1$thermo_lab <- as.numeric(mapvalues(us1$thermo_lab, c("Extremely cold0","1","2","3","4",
                                                         "Neither warm nor cold5","6","7","8",
                                                         "9","Extremely warm10"), 0:10))
us1$thermo_lib <- as.numeric(mapvalues(us1$thermo_lib, c("Extremely cold0","1","2","3","4",
                                                         "Neither warm nor cold5","6","7","8",
                                                         "9","Extremely warm10"), 0:10))
us1$thermo_diff <- us1$thermo_lib - us1$thermo_lab



# Produce binary measure of racial identity 

us1$White <- ifelse(us1$Race_Class_Identity == "White", 1, 0)
us1$White = factor(us1$White, levels = c(0, 1), labels = c('Minority', 'White'))


## wave 3

# Clean redistribution variable - change item names, recode to numeric, produce mean

us3$redist1 <- us3$Q4_1
us3$redist2 <- us3$Q4_2
us3$redist3 <- us3$Q4_3
us3$redist4 <- us3$Q4_4
us3$redist5 <- us3$Q4_5
us3$redist6 <- us3$Q4_6
us3$redist1 <- as.numeric(mapvalues(us3$redist1,c("Strongly disagree","Somewhat disagree","Neither agree nor disagree",
                                                  "Somewhat agree","Strongly agree"),
                                    1:5))
us3$redist2 <- as.numeric(mapvalues(us3$redist2,c("Strongly disagree","Somewhat disagree","Neither agree nor disagree",
                                                  "Somewhat agree","Strongly agree"),
                                    1:5))
us3$redist3 <- as.numeric(mapvalues(us3$redist3,c("Strongly disagree","Somewhat disagree","Neither agree nor disagree",
                                                  "Somewhat agree","Strongly agree"),
                                    1:5))
us3$redist4 <- as.numeric(mapvalues(us3$redist4,c("Strongly disagree","Somewhat disagree","Neither agree nor disagree",
                                                  "Somewhat agree","Strongly agree"),
                                    5:1))
us3$redist5 <- as.numeric(mapvalues(us3$redist5,c("Strongly disagree","Somewhat disagree","Neither agree nor disagree",
                                                  "Somewhat agree","Strongly agree"),
                                    1:5))
us3$redist6 <- as.numeric(mapvalues(us3$redist6,c("Strongly disagree","Somewhat disagree","Neither agree nor disagree",
                                                  "Somewhat agree","Strongly agree"),
                                    1:5))
us3$redist_avg <- apply(us3[,paste("redist",1:6,sep="")],1,mean,na.rm=T)

us3$redist_fac <- fa(us3[,paste("redist",1:6,sep="")],nfactors=1,rotate="Promax", cor = "poly")$scores[,1]
us3$redist_fac_1to5 <- scales:::rescale(us3$redist_fac,to=c(1,5))


# Clean elite domination variable - change item names, recode to numeric, produce mean and factor score

us3$rigged1 <- us3$Q19_1
us3$rigged2 <- us3$Q19_2
us3$rigged3 <- us3$Q19_3
us3$rigged4 <- us3$Q19_4
us3$rigged5 <- us3$Q19_5
us3$rigged6 <- us3$Q19_6
us3$rigged1 <- as.numeric(mapvalues(us3$rigged1,c("Strongly agree","Somewhat agree",
                                                  "Neither agree nor disagree","Somewhat disagree",
                                                  "Strongly disagree"), 5:1))
us3$rigged2 <- as.numeric(mapvalues(us3$rigged2,c("Strongly agree","Somewhat agree",
                                                  "Neither agree nor disagree","Somewhat disagree",
                                                  "Strongly disagree"), 5:1))
us3$rigged3 <- as.numeric(mapvalues(us3$rigged3,c("Strongly agree","Somewhat agree",
                                                  "Neither agree nor disagree","Somewhat disagree",
                                                  "Strongly disagree"), 5:1))
us3$rigged4 <- as.numeric(mapvalues(us3$rigged4,c("Strongly agree","Somewhat agree",
                                                  "Neither agree nor disagree","Somewhat disagree",
                                                  "Strongly disagree"), 5:1))
us3$rigged5 <- as.numeric(mapvalues(us3$rigged5,c("Strongly agree","Somewhat agree",
                                                  "Neither agree nor disagree","Somewhat disagree",
                                                  "Strongly disagree"), 1:5))
us3$rigged6 <- as.numeric(mapvalues(us3$rigged6,c("Strongly agree","Somewhat agree",
                                                  "Neither agree nor disagree","Somewhat disagree",
                                                  "Strongly disagree"), 5:1))

us3$elitedom <- apply(us3[,paste("rigged",1:6,sep="")],1,mean,na.rm=T)

us3$elitedom_polychoric_factor_score <- fa(us3[,paste("rigged",1:6,sep="")],nfactors=1,rotate="Promax", cor = "poly")$scores[,1]
us3$elitedom_polychoric_factor_score <- scales:::rescale(us3$elitedom_polychoric_factor_score,to=c(1,5))








#######################
#### Organize variables: Switzerland
#######################

## wave 1

# Create average measure of redistribution attitudes

ch1$redist_avg_w1 <- apply(ch1[,paste("redist",1:6,sep="")],1,mean,na.rm=T)



# Form binary variable of news consumption

ch1$news_frequency <- ifelse(ch1$news_freq_tel == "6" | 
                               ch1$news_freq_tel == "7" |
                               ch1$news_freq_paper == "6" | 
                               ch1$news_freq_paper == "7" |
                               ch1$news_freq_radio == "6" | 
                               ch1$news_freq_radio == "7" |
                               ch1$news_freq_social == "6" | 
                               ch1$news_freq_social == "7" |
                               ch1$news_freq_other == "6" | 
                               ch1$news_freq_other == "7", 1, 0)

ch1$news_frequency[is.na(ch1$news_frequency)] <- 0
ch1$news_frequency = factor(ch1$news_frequency, levels = c(0, 1), 
                            labels = c('Low Exposure to Financial News', 'High Exposure to Financial News'))



# Produce mean and factor score of racial resentment variable

ch1$rr_avg <- apply(ch1[,paste("rr",1:4,sep="")],1,mean,na.rm=T)


# Produce mean and factor score of individualism scale items

ch1$indiv_avg <- apply(ch1[,paste("indiv",1:5,sep="")],1,mean,na.rm=T)

ch1$indiv_avg_polychoric_factor <- fa(ch1[,paste(c("indiv1","indiv2", "indiv3", "indiv4", "indiv5"))],nfactors=1,rotate="Promax", cor = "poly")$scores[,1]
ch1$indiv_avg_polychoric_factor <- scales:::rescale(ch1$indiv_avg_polychoric_factor,to=c(1,5))


# Produce affective polarization / partisan identification construct 

ch1$thermo_lab <- ch1$Q30a_g_1 # SP
ch1$thermo_lib <- ch1$Q30a_g_2 # SVP
ch1$thermo_lab <- as.numeric(mapvalues(ch1$thermo_lab, c("0","1","2","3","4",
                                                         "5","6","7","8",
                                                         "9","10"), 0:10))
ch1$thermo_lib <- as.numeric(mapvalues(ch1$thermo_lib, c("0","1","2","3","4",
                                                         "5","6","7","8",
                                                         "9","10"), 0:10))
ch1$thermo_diff <- ch1$thermo_lib - ch1$thermo_lab


# Add factor score

ch1$rr1 <- (6-ch1$rr1)
ch1$rr2 <- (6-ch1$rr2)
ch1$rr3 <- (6-ch1$rr3)
ch1$rr4 <- (6-ch1$rr4)

ch1$rr_avg_polychoric_factor <- fa(ch1[,paste(c("rr1","rr2", "rr3", "rr4"))],nfactors=1,rotate="Promax", cor = "poly")$scores[,1]
ch1$rr_avg_polychoric_factor <- scales:::rescale(ch1$rr_avg_polychoric_factor,to=c(1,5))

## wave 3

# Produce mean of redistribution variable

ch3$redist_avg <- apply(ch3[,paste("redist",1:6,sep="")],1,mean,na.rm=T)

ch3$redist_fac <- fa(ch3[,paste("redist",1:6,sep="")],nfactors=1,rotate="Promax", cor = "poly")$scores[,1]
ch3$redist_fac_1to5 <- scales:::rescale(ch3$redist_fac,to=c(1,5))


# Produce mean and factor score of elite domination variable


ch3$elitedom <- apply(ch3[,paste("rigged",1:6,sep="")],1,mean,na.rm=T)

ch3$elitedom_polychoric_factor_score <- fa(ch3[,paste("rigged",1:6,sep="")],nfactors=1,rotate="Promax", cor = "poly")$scores[,1]
ch3$elitedom_polychoric_factor_score <- scales:::rescale(ch3$elitedom_polychoric_factor_score,to=c(1,5))





############
############
#Combining Waves 1 and 3 for each country, and adding a country identifier variable
############
############

au_combined <-  merge(au1, au3[,c("redist_avg", "elitedom", "elitedom_polychoric_factor_score", "redist_fac_1to5",
                                  "treatment", "longitudinal_ID", "weight3",
                                  "rigged1", "rigged2", "rigged3", "rigged4", "rigged5", "rigged6")],
                      by="longitudinal_ID",all=T)
au_combined$country <- "AU"

fr_combined <-  merge(fr1, fr3[,c("redist_avg", "elitedom", "elitedom_polychoric_factor_score", "redist_fac_1to5",
                                  "treatment", "longitudinal_ID", "weight3",
                                  "rigged1", "rigged2", "rigged3", "rigged4", "rigged5", "rigged6")],
                      by="longitudinal_ID",all=T)
fr_combined$country <- "FR"

de_combined <-  merge(de1, de3[,c("redist_avg", "elitedom", "elitedom_polychoric_factor_score", "redist_fac_1to5",
                                  "treatment", "longitudinal_ID", "weight3",
                                  "rigged1", "rigged2", "rigged3", "rigged4", "rigged5", "rigged6")],
                      by="longitudinal_ID",all=T)
de_combined$country <- "DE"

uk_combined <-  merge(uk1, uk3[,c("redist_avg", "elitedom", "elitedom_polychoric_factor_score", "redist_fac_1to5",
                                  "treatment", "longitudinal_ID", "weight3",
                                  "rigged1", "rigged2", "rigged3", "rigged4", "rigged5", "rigged6")],
                      by="longitudinal_ID",all=T)
uk_combined$country <- "UK"

us_combined <-  merge(us1, us3[,c("redist_avg", "elitedom", "elitedom_polychoric_factor_score", "redist_fac_1to5",
                                  "treatment", "longitudinal_ID", "weight3",
                                  "rigged1", "rigged2", "rigged3", "rigged4", "rigged5", "rigged6")],
                      by="longitudinal_ID",all=T)
us_combined$country <- "US"

ch1$Q30a_g_1 <- as.character(ch1$Q30a_g_1) #fixing variable type to main consistency with other countries
ch1$Q30a_g_2 <- as.character(ch1$Q30a_g_2)

ch_combined <-  merge(ch1, ch3[,c("redist_avg", "elitedom", "elitedom_polychoric_factor_score", "redist_fac_1to5",
                                  "treatment", "longitudinal_ID", "weight3",
                                  "rigged1", "rigged2", "rigged3", "rigged4", "rigged5", "rigged6")],
                      by="longitudinal_ID",all=T)
ch_combined$country <- "CH"


############
# Subsetting into Treatment Groups
############

au_nb <- subset(au_combined, treatment=="Non-Bank Article")
au_rig <- subset(au_combined, treatment=="Rigged Economy Article")
ch_nb <- subset(ch_combined, treatment=="Non-Bank Article")
ch_rig <- subset(ch_combined, treatment=="Rigged Economy Article")
de_nb <- subset(de_combined, treatment=="Non-Bank Article")
de_rig <- subset(de_combined, treatment=="Rigged Economy Article")
fr_nb <- subset(fr_combined, treatment=="Non-Bank Article")
fr_rig <- subset(fr_combined, treatment=="Rigged Economy Article")
uk_nb <- subset(uk_combined, treatment=="Non-Bank Article")
uk_rig <- subset(uk_combined, treatment=="Rigged Economy Article")
us_nb <- subset(us_combined, treatment=="Non-Bank Article")
us_rig <- subset(us_combined, treatment=="Rigged Economy Article")


############
############
# Combining All Countries and Waves into Single Dataset
############
############

# minor edits to make the variables consistent and permit bind_rows to function

ch_combined$female = factor(ch_combined$female, levels = c(0, 1), labels = c('Male', 'Female'))
uk_combined$female = factor(uk_combined$female, levels = c(0, 1), labels = c('Male', 'Female'))



# Combine all countries 

combo = bind_rows(au_combined, us_combined, uk_combined, de_combined, fr_combined, ch_combined)

combo <- subset(combo, combo$treatment == "Non-Bank Article" | combo$treatment == "Rigged Economy Article")





#######################
#######################
#### Figure 3a     #### 
#######################
#######################

# Store coefficients and standard errors from the analytical models:

diff_au_rigged <- -wtd.t.test(au_nb$redist_avg,au_rig$redist_avg,weight=au_nb$weight3,weighty=au_rig$weight3,samedata=F)$additional[1]
diff_au_rigged_se <- wtd.t.test(au_nb$redist_avg,au_rig$redist_avg,weight=au_nb$weight3,weighty=au_rig$weight3,samedata=F)$additional[4]

diff_fr_rigged <- -wtd.t.test(fr_nb$redist_avg,fr_rig$redist_avg,weight=fr_nb$weight3,weighty=fr_rig$weight3,samedata=F)$additional[1]
diff_fr_rigged_se <- wtd.t.test(fr_nb$redist_avg,fr_rig$redist_avg,weight=fr_nb$weight3,weighty=fr_rig$weight3,samedata=F)$additional[4]

diff_de_rigged <- -wtd.t.test(de_nb$redist_avg,de_rig$redist_avg,weight=de_nb$weight3,weighty=de_rig$weight3,samedata=F)$additional[1]
diff_de_rigged_se <- wtd.t.test(de_nb$redist_avg,de_rig$redist_avg,weight=de_nb$weight3,weighty=de_rig$weight3,samedata=F)$additional[4]

diff_ch_rigged <- -wtd.t.test(ch_nb$redist_avg,ch_rig$redist_avg,weight=ch_nb$weight3,weighty=ch_rig$weight3,samedata=F)$additional[1]
diff_ch_rigged_se <- wtd.t.test(ch_nb$redist_avg,ch_rig$redist_avg,weight=ch_nb$weight3,weighty=ch_rig$weight3,samedata=F)$additional[4]

diff_uk_rigged <- -wtd.t.test(uk_nb$redist_avg,uk_rig$redist_avg,weight=uk_nb$weight3,weighty=uk_rig$weight3,samedata=F)$additional[1]
diff_uk_rigged_se <- wtd.t.test(uk_nb$redist_avg,uk_rig$redist_avg,weight=uk_nb$weight3,weighty=uk_rig$weight3,samedata=F)$additional[4]

diff_us_rigged <- -wtd.t.test(us_nb$redist_avg,us_rig$redist_avg,weight=us_nb$weight3,weighty=us_rig$weight3,samedata=F)$additional[1]
diff_us_rigged_se <- wtd.t.test(us_nb$redist_avg,us_rig$redist_avg,weight=us_nb$weight3,weighty=us_rig$weight3,samedata=F)$additional[4]


# Plot

par(mar=c(6.1,5.1,2,2.1))
plot(x=c(),y=c(),ylim=c(-0.3,0.3),xlim=c(0.1,0.8), ylab="Treatment Effect", main="", xaxt="n", xlab="",cex.lab=1.4)
points(0.2,diff_au_rigged,pch=16)
lines(x=c(0.2,0.2),y=c(diff_au_rigged-1.64*diff_au_rigged_se,diff_au_rigged+1.64*diff_au_rigged_se),lwd=2)
points(0.3,diff_fr_rigged,pch=16)
lines(x=c(0.3,0.3),y=c(diff_fr_rigged-1.64*diff_fr_rigged_se,diff_fr_rigged+1.64*diff_fr_rigged_se),lwd=2)
points(0.4,diff_de_rigged,pch=16)
lines(x=c(0.4,0.4),y=c(diff_de_rigged-1.64*diff_de_rigged_se,diff_de_rigged+1.64*diff_de_rigged_se),lwd=2)
points(0.5,diff_ch_rigged,pch=16)
lines(x=c(0.5,0.5),y=c(diff_ch_rigged-1.64*diff_ch_rigged_se,diff_ch_rigged+1.64*diff_ch_rigged_se),lwd=2)
points(0.6,diff_uk_rigged,pch=16)
lines(x=c(0.6,0.6),y=c(diff_uk_rigged-1.64*diff_uk_rigged_se,diff_uk_rigged+1.64*diff_uk_rigged_se),lwd=2)
points(0.7,diff_us_rigged,pch=16)
lines(x=c(0.7,0.7),y=c(diff_us_rigged-1.64*diff_us_rigged_se,diff_us_rigged+1.64*diff_us_rigged_se),lwd=2)
abline(h=0,lty="dashed")
axis(1, at=seq(0.2,0.7,by=0.1), labels = FALSE)
text(seq(0.03,0.73,by=0.1), par("usr")[3] - 0.10, labels = c("","Australia","France","Germany","Switzerland", "United Kingdom","United States",""), srt = 45, pos = 1, xpd = TRUE,cex=1.2)


# Effect Sizes as Noted in Manuscript

diff_au_rigged #Australian Effect Size
diff_ch_rigged #Swiss Effect Size
diff_de_rigged #German Effect Size
diff_fr_rigged #French Effect Size
diff_uk_rigged #UK Effect Size
diff_us_rigged #US Effect Size




#######################
#######################
#### Figure 3b     #### 
#######################
#######################

#### Generate Functions

weighted.average <- function(x, w){
  ## Sum of the weights 
  sum.w <- sum(w, na.rm = T)
  ## Sum of the weighted $x_i$ 
  xw <- sum(w*x, na.rm = T)
  
  ## Return the weighted average 
  return(xw/sum.w)
}

weighted.se.mean <- function(x, w, na.rm = T){
  ## Remove NAs 
  if (na.rm) {
    i <- !is.na(x)
    w <- w[i]
    x <- x[i]
  }
  
  ## Calculate effective N and correction factor
  n_eff <- (sum(w))^2/(sum(w^2))
  correction = n_eff/(n_eff-1)
  
  ## Get weighted variance 
  numerator = sum(w*(x-wtd.mean(x,w))^2)
  denominator = sum(w)
  
  ## get weighted standard error of the mean 
  se_x = sqrt((correction * (numerator/denominator))/n_eff)
  return(se_x)
}

summaries <- data.frame(
  #Country=c("AU","AU","FR","FR","DE","DE","CH","CH","UK","UK","US","US"),
  Country=c("Australia","Australia","France","France","Germany","Germany","Switzerland","Switzerland","United Kingdom","United Kingdom","United States","United States"),
  Condition=rep(c("Control","Treatment"),6),
  wtd.mean=c(wtd.mean(au_nb$redist_avg,au_nb$weight3),wtd.mean(au_rig$redist_avg,au_rig$weight3),
             wtd.mean(fr_nb$redist_avg,fr_nb$weight3),wtd.mean(fr_rig$redist_avg,fr_rig$weight3),
             wtd.mean(de_nb$redist_avg,de_nb$weight3),wtd.mean(de_rig$redist_avg,de_rig$weight3),
             wtd.mean(ch_nb$redist_avg,ch_nb$weight3),wtd.mean(ch_rig$redist_avg,ch_rig$weight3),
             wtd.mean(uk_nb$redist_avg,uk_nb$weight3),wtd.mean(uk_rig$redist_avg,uk_rig$weight3),
             wtd.mean(us_nb$redist_avg,us_nb$weight3),wtd.mean(us_rig$redist_avg,us_rig$weight3)),
  wtd.se=c(weighted.se.mean(au_nb$redist_avg,au_nb$weight3),weighted.se.mean(au_rig$redist_avg,au_rig$weight3),
           weighted.se.mean(fr_nb$redist_avg,fr_nb$weight3),weighted.se.mean(fr_rig$redist_avg,fr_rig$weight3),
           weighted.se.mean(de_nb$redist_avg,de_nb$weight3),weighted.se.mean(de_rig$redist_avg,de_rig$weight3),
           weighted.se.mean(ch_nb$redist_avg,ch_nb$weight3),weighted.se.mean(ch_rig$redist_avg,ch_rig$weight3),
           weighted.se.mean(uk_nb$redist_avg,uk_nb$weight3),weighted.se.mean(uk_rig$redist_avg,uk_rig$weight3),
           weighted.se.mean(us_nb$redist_avg,us_nb$weight3),weighted.se.mean(us_rig$redist_avg,us_rig$weight3))
)




#### Plot

ggplot(summaries, aes(x = factor(Country), y = wtd.mean, fill = Condition), ylim=c(1,5)) +
  geom_bar(stat = "identity", position = "dodge") +
  geom_errorbar(aes(ymax = wtd.mean + wtd.se, ymin = wtd.mean - wtd.se),
                position = position_dodge(0.9), width = 0.25, color = "Gray25") +
  ylab("Weighted Mean\n") + 
  theme(axis.title.x=element_blank(),panel.background = element_rect(fill = "white", colour = "black"),
        axis.text=element_text(size=14),
        axis.text.x=element_text(angle=45, vjust=1, hjust=1.1, color="black"),
        axis.text.y=element_text(color="black", size = 12),
        axis.title=element_text(size=14),
        legend.title = element_text(size=12),
        legend.text = element_text(size=12)) +
  scale_fill_grey(start=0.6) + 
  coord_cartesian(ylim = c(3, 4))





#######################
#######################
#### Figure 4  #### 
#######################
#######################

# Add variable that distinguishes between the US and the five-country bloc

combo$FiveCountries <- ifelse(combo$country == "AU" | combo$country == "FR" | combo$country == "DE" | combo$country == "UK" | combo$country == "CH", "Yes", "No")



######## Five Country Bloc ########

model.m <- lm(elitedom_polychoric_factor_score ~ treatment + age + female + income + educ_fourlevels + lr + country, data = subset(combo, FiveCountries == "Yes"))
model.y <- lm(redist_avg ~ elitedom_polychoric_factor_score + treatment + age + female + income + educ_fourlevels + lr + country, data = subset(combo, FiveCountries == "Yes"))

out.1 <- mediate(model.m, model.y, sims = 1000, boot = TRUE, treat = "treatment", mediator = "elitedom_polychoric_factor_score")


summary(out.1)
plot(out.1)
summary(model.m)
summary(model.y)

### Sensitivity Analysis -  ###

med.cont1 <- mediate(model.m, model.y, sims=1000, treat = "treatment", mediator = "elitedom_polychoric_factor_score")
sens.cont1 <- medsens(med.cont1, rho.by = 0.01)

summary(sens.cont1)

plot(sens.cont1, sens.par = "rho")
plot(sens.cont1, sens.par = "rho", xlim = c(-.9,.9), ylim = c(-.9,.9))

### Combined Plot - Mediation and Sensitivity Analyses -  ###

par(mfrow=c(1,2))
plot(out.1,main="Mediation Analysis: \nFive-Country Bloc")
plot(sens.cont1,main="Sensitivity Analysis: \nFive-Country Bloc")
legend("topright",legend= paste("rho:",sens.cont1$err.cr.d),bty="n")




######## United States ########

model.m <- lm(elitedom_polychoric_factor_score ~ treatment + age + female + income + educ_fourlevels + lr, data = subset(combo, country == "US"))
model.y <- lm(redist_avg ~ elitedom_polychoric_factor_score + treatment + age + female + income + educ_fourlevels + lr, data = subset(combo, country == "US"))

out.1 <- mediate(model.m, model.y, sims = 1000, boot = TRUE, treat = "treatment", mediator = "elitedom_polychoric_factor_score")


summary(out.1)
plot(out.1)
summary(model.m)
summary(model.y)

### Sensitivity Analysis -  ###

med.cont1 <- mediate(model.m, model.y, sims=1000, treat = "treatment", mediator = "elitedom_polychoric_factor_score")
sens.cont1 <- medsens(med.cont1, rho.by = 0.01)

summary(sens.cont1)

plot(sens.cont1, sens.par = "rho")
plot(sens.cont1, sens.par = "rho", xlim = c(-.9,.9), ylim = c(-.9,.9))

### Combined Plot - Mediation and Sensitivity Analyses -  ###

par(mfrow=c(1,2))
plot(out.1,main="Mediation Analysis: \nUnited States")
plot(sens.cont1,main="Sensitivity Analysis: \nUnited States")
legend("topright",legend= paste("rho:",sens.cont1$err.cr.d),bty="n")
par(mfrow=c(1,1))




#######################
#######################
#### Figure 5a     #### 
#######################
#######################

# Modeling interaction of treatment condition and political orientation
mod2 <- lm(redist_avg ~ treatment*lr, data=subset(combo, country=="US"), weights = weight3)
summary(mod2)

# Plotting interaction effeft
par(mfrow=c(1,1),oma=c(1,1,0,0),mai=c(0.5,0.4,0,0),mar=c(4,4,1,1)) 
thermo_sim <- 0:10

treat <- coef(mod2)[2] + coef(mod2)[4] * thermo_sim
treat_se <- sqrt(vcov(mod2)[2, 2] + thermo_sim^2 * vcov(mod2)[4, 4] + 2 * thermo_sim * vcov(mod2)[2, 4])
plot(x = thermo_sim, y = treat, type = "p",pch=16,
     xlab = "Political Orientation",
     ylab = "Effect of Riged System Treatment",
     ylim = c(-0.4,0.4), cex.lab=1.1)
for(i in 1:11){
  lines(c(i-1,i-1),c(treat[i],treat[i] + 1.96 * treat_se[i]),lwd=2)
  lines(c(i-1,i-1),c(treat[i],treat[i] - 1.96 * treat_se[i]),lwd=2)
}
abline(h=0,lty="dashed")

par(new=T)
hist(subset(combo, country == "US")$lr,xlab=NULL,ylab=NULL,axes=F,main=NULL,freq=F,col=rgb(0.1,0.1,0.1,0.1),ylim=c(0,2),breaks = c(-1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10))




#######################
#######################
#### Figure 5b     #### 
#######################
#######################

# Modeling interaction of treatment condition and political orientation
mod2 <- lm(redist_avg ~ treatment*lr, data=subset(combo, FiveCountries=="Yes"), weights = weight3)
summary(mod2)

# Plotting interaction effeft
par(mfrow=c(1,1),oma=c(1,1,0,0),mai=c(0.5,0.4,0,0),mar=c(4,4,1,1)) 
thermo_sim <- 0:10

treat <- coef(mod2)[2] + coef(mod2)[4] * thermo_sim
treat_se <- sqrt(vcov(mod2)[2, 2] + thermo_sim^2 * vcov(mod2)[4, 4] + 2 * thermo_sim * vcov(mod2)[2, 4])
plot(x = thermo_sim, y = treat, type = "p",pch=16,
     xlab = "Political Orientation",
     ylab = "Effect of Riged System Treatment",
     ylim = c(-0.4,0.4), cex.lab=1.1)
for(i in 1:11){
  lines(c(i-1,i-1),c(treat[i],treat[i] + 1.96 * treat_se[i]),lwd=2)
  lines(c(i-1,i-1),c(treat[i],treat[i] - 1.96 * treat_se[i]),lwd=2)
}
abline(h=0,lty="dashed")

par(new=T)
hist(subset(combo, FiveCountries == "Yes")$lr,xlab=NULL,ylab=NULL,axes=F,main=NULL,freq=F,col=rgb(0.1,0.1,0.1,0.1),ylim=c(0,2),breaks = c(-1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10))






#######################
#######################
#### Figure 6      #### 
#######################
#######################

# Normalizing Variables to ease the comparison of effect sizes #

combo_stand <- combo
combo_stand$lr_stand <- rescale(combo_stand$lr)
combo_stand$income_stand <- rescale(combo_stand$income)
combo_stand$elitedom_polychoric_factor_score_stand <- rescale(combo_stand$elitedom_polychoric_factor_score)
combo_stand$trust_govt_w1_stand <- rescale(combo_stand$trust_govt_w1)
combo_stand$rr_avg_polychoric_factor_stand <- rescale(combo_stand$rr_avg_polychoric_factor)
combo_stand$indiv_avg_polychoric_factor_stand <- rescale(combo_stand$indiv_avg_polychoric_factor)
combo_stand$lf4_stand <- rescale(combo_stand$lf4)
combo_stand$efficacy_avg_stand <- rescale(combo_stand$efficacy_avg)
combo_stand$educ_fourlevels = as.factor(combo_stand$educ_fourlevels)

# Converting age variable to ordinal age_groups variable, and fixing reversed variable

combo$age_groups <- ifelse(combo$age < 35, 0,
                           ifelse(combo$age > 34 & combo$age < 50, 1,
                                  ifelse(combo$age>49 & combo$age < 63, 2,3)))
combo$age_groups = factor(combo$age_groups, levels = c(0, 1, 2, 3), labels = c('Below 35', '35 - 50', '50 - 62', '63+'))
combo_stand$age_groups <- combo$age_groups

###########################
# Generating Country Level Regressions #
###########################

PredictingAU_stand <- lm(redist_avg_w1 ~
                           lr_stand + news_frequency + female + age_groups + educ_fourlevels + income_stand
                         + trust_govt_w1_stand + rr_avg_polychoric_factor_stand + indiv_avg_polychoric_factor_stand
                         + lf4_stand  + efficacy_avg_stand
                         , data = subset(combo_stand, country=="AU"))

PredictingCH_stand <- lm(redist_avg_w1 ~
                           lr_stand + news_frequency + female + age_groups + educ_fourlevels + income_stand
                         + trust_govt_w1_stand + rr_avg_polychoric_factor_stand + indiv_avg_polychoric_factor_stand
                         + lf4_stand  + efficacy_avg_stand
                         , data = subset(combo_stand, country=="CH"))
summary(PredictingCH_stand)

PredictingDE_stand <- lm(redist_avg_w1 ~
                           lr_stand + news_frequency + female + age_groups + educ_fourlevels + income_stand
                         + trust_govt_w1_stand + rr_avg_polychoric_factor_stand + indiv_avg_polychoric_factor_stand
                         + lf4_stand  + efficacy_avg_stand
                         , data = subset(combo_stand, country=="DE"))

PredictingFR_stand <- lm(redist_avg_w1 ~
                           lr_stand + news_frequency + female + age_groups + educ_fourlevels + income_stand
                         + trust_govt_w1_stand + rr_avg_polychoric_factor_stand + indiv_avg_polychoric_factor_stand
                         + lf4_stand  + efficacy_avg_stand
                         , data = subset(combo_stand, country=="FR"))

PredictingUK_stand <- lm(redist_avg_w1 ~
                           lr_stand + news_frequency + female + age_groups + educ_fourlevels + income_stand
                         + trust_govt_w1_stand + rr_avg_polychoric_factor_stand + indiv_avg_polychoric_factor_stand
                         + lf4_stand  + efficacy_avg_stand
                         , data = subset(combo_stand, country=="UK"))

PredictingUS_stand <- lm(redist_avg_w1 ~
                           lr_stand + news_frequency + female + age_groups + educ_fourlevels + income_stand
                         + trust_govt_w1_stand + rr_avg_polychoric_factor_stand + indiv_avg_polychoric_factor_stand
                         + lf4_stand  + efficacy_avg_stand
                         , data = subset(combo_stand, country=="US"))



###########################
# Plotting Figure 6
###########################

plot_summs(PredictingAU_stand, PredictingCH_stand, PredictingDE_stand, PredictingFR_stand, PredictingUK_stand, PredictingUS_stand,
           robust = TRUE,
           legend.title = "Country",
           model.names = c("AU", "CH", "DE", "FR", "UK", "US"),
           coefs = c("Political Orientation" = "lr_stand",
                     "News Consumption" = "news_frequencyHigh Exposure to Financial News",
                     "Gender (Female)" = "femaleFemale",
                     "Age Group 35-50" = "age_groups35 - 50",
                     "Age Group 50-62" = "age_groups50 - 62",
                     "Age Group 63+" = "age_groups63+",
                     "Education Level 2" = "educ_fourlevels2",
                     "Education Level 3" = "educ_fourlevels3",
                     "Education Level 4" = "educ_fourlevels4",
                     "Income" = "income_stand",
                     "Trust in Government" = "trust_govt_w1_stand",
                     "Racial Resentment" = "rr_avg_polychoric_factor_stand",
                     "Economic Individualism" = "indiv_avg_polychoric_factor_stand",
                     "Government Inefficiency" = "lf4_stand",
                     "Political Efficacy" = "efficacy_avg_stand"))







#######################
#######################
#### Figure 7  #### 
#######################
#######################

######## United States ########

# Modeling Effect

mod2 <- lm(redist_avg ~ treatment*lf4, data=subset(combo, country=="US"), weights = weight3) 
summary(mod2)

# Plotting Effect

par(mfrow=c(1,1),oma=c(1,1,0,0),mai=c(0.5,0.4,0,0),mar=c(4,4,1,1)) # oma is originally 0 0 0 0. mai is 0.6732 0.5412 0.5412 0.2772
thermo_sim <- 1:5

treat <- coef(mod2)[2] + coef(mod2)[4] * thermo_sim
treat_se <- sqrt(vcov(mod2)[2, 2] + thermo_sim^2 * vcov(mod2)[4, 4] + 2 * thermo_sim * vcov(mod2)[2, 4])
plot(x = thermo_sim, y = treat, type = "p",pch=16,
     xlab = "Government Intervention Leads to Too Much Red Tape",
     ylab = "Treatment Effect",
     ylim = c(-1.0,1.2), cex.lab=1.1)
for(i in 1:5){
  lines(c(i,i),c(treat[i],treat[i] + 1.96 * treat_se[i]),lwd=2)
  lines(c(i,i),c(treat[i],treat[i] - 1.96 * treat_se[i]),lwd=2)
}
abline(h=0,lty="dashed")

par(new=T)
hist(subset(combo, country == "US")$lf4,xlab=NULL,ylab=NULL,axes=F,main=NULL,freq=F,col=rgb(0.1,0.1,0.1,0.1),ylim=c(0,2.5), breaks = c(0, 1, 2, 3, 4, 5))


######## Five-Country Bloc ########

# Modeling Effect

mod2 <- lm(redist_avg ~ treatment*lf4, data=subset(combo, FiveCountries == "Yes"), weights = weight3)       
summary(mod2)


# Plotting Effect

par(mfrow=c(1,1),oma=c(1,1,0,0),mai=c(0.5,0.4,0,0),mar=c(4,4,1,1)) # oma is originally 0 0 0 0. mai is 0.6732 0.5412 0.5412 0.2772
thermo_sim <- 1:5

treat <- coef(mod2)[2] + coef(mod2)[4] * thermo_sim
treat_se <- sqrt(vcov(mod2)[2, 2] + thermo_sim^2 * vcov(mod2)[4, 4] + 2 * thermo_sim * vcov(mod2)[2, 4])
plot(x = thermo_sim, y = treat, type = "p",pch=16,
     xlab = "Government Intervention Leads to Too Much Red Tape",
     ylab = "Treatment Effect",
     ylim = c(-1.0,1.2), cex.lab=1.1)
for(i in 1:5){
  lines(c(i,i),c(treat[i],treat[i] + 1.96 * treat_se[i]),lwd=2)
  lines(c(i,i),c(treat[i],treat[i] - 1.96 * treat_se[i]),lwd=2)
}
abline(h=0,lty="dashed")

par(new=T)
hist(subset(combo, FiveCountries == "Yes")$lf4,xlab=NULL,ylab=NULL,axes=F,main=NULL,freq=F,col=rgb(0.1,0.1,0.1,0.1),ylim=c(0,2.5), breaks = c(0, 1, 2, 3, 4, 5))


######## Analysis of Three-Way Interaction Effect ########

summary(lm(redist_avg ~ treatment*lf4*FiveCountries, data=combo, weights = weight3))





###################################
###################################
####### Online Appendix A  ########
###################################
###################################

# Function Code for Balance Test

get_mean_sd <- function(vector){
  mean <- mean(vector, na.rm = T)
  sd <- sd(vector, na.rm = T)
  sem <- sd/(sqrt(length(vector)))
  output <- paste(round(mean, 3), " (", round(sem, 3), ")", collapse = "", sep = "")
  return(output)
}


### Table A1

balance_data <- group_by(combo, treatment) %>%
  dplyr::summarize(Age = get_mean_sd(age),
            Gender = get_mean_sd(as.numeric(female)), 
            Political_Orientation = get_mean_sd(as.numeric(lr)),
            Education = get_mean_sd(as.numeric(educ_fourlevels)),
            Income = get_mean_sd(as.numeric(income)))
balance_data



### Table A2

# Convert treatment group into numeric variable

combo$treatment_numeric <- ifelse(combo$treatment == "Rigged Economy Article", 1, 0)

# Run regression

ProbitBalanceCheck <- glm(treatment_numeric ~  
                              age + as.numeric(female) + lr + as.numeric(educ_fourlevels) + income, data=combo) 
summary(ProbitBalanceCheck)
logLik(ProbitBalanceCheck)




###################################
###################################
####### Online Appendix C  ########
###################################
###################################

# Factor Loadings

fa(subset(au3, treatment == "Rigged Economy Article" | treatment == "Non-Bank Article")
          [,paste("rigged",1:6,sep="")],nfactors=1,rotate="Promax", cor = "poly")

fa(subset(de3, treatment == "Rigged Economy Article" | treatment == "Non-Bank Article")
   [,paste("rigged",1:6,sep="")],nfactors=1,rotate="Promax", cor = "poly")

fa(subset(fr3, treatment == "Rigged Economy Article" | treatment == "Non-Bank Article")
   [,paste("rigged",1:6,sep="")],nfactors=1,rotate="Promax", cor = "poly")

fa(subset(uk3, treatment == "Rigged Economy Article" | treatment == "Non-Bank Article")
   [,paste("rigged",1:6,sep="")],nfactors=1,rotate="Promax", cor = "poly")

fa(subset(us3, treatment == "Rigged Economy Article" | treatment == "Non-Bank Article")
   [,paste("rigged",1:6,sep="")],nfactors=1,rotate="Promax", cor = "poly")

fa(subset(ch3, treatment == "Rigged Economy Article" | treatment == "Non-Bank Article")
   [,paste("rigged",1:6,sep="")],nfactors=1,rotate="Promax", cor = "poly")

fa(subset(combo, treatment == "Rigged Economy Article" | treatment == "Non-Bank Article")
   [,paste("rigged",1:6,sep="")],nfactors=1,rotate="Promax", cor = "poly")


# Cronbach's Alpha Scores

summary(psych::alpha(subset(au3, treatment == "Rigged Economy Article" | treatment == "Non-Bank Article")[,paste("rigged",1:6,sep="")]))
summary(psych::alpha(subset(ch3, treatment == "Rigged Economy Article" | treatment == "Non-Bank Article")[,paste("rigged",1:6,sep="")]))
summary(psych::alpha(subset(de3, treatment == "Rigged Economy Article" | treatment == "Non-Bank Article")[,paste("rigged",1:6,sep="")]))
summary(psych::alpha(subset(fr3, treatment == "Rigged Economy Article" | treatment == "Non-Bank Article")[,paste("rigged",1:6,sep="")]))
summary(psych::alpha(subset(uk3, treatment == "Rigged Economy Article" | treatment == "Non-Bank Article")[,paste("rigged",1:6,sep="")]))
summary(psych::alpha(subset(us3, treatment == "Rigged Economy Article" | treatment == "Non-Bank Article")[,paste("rigged",1:6,sep="")]))

# Inter-Item Correlations

interitem_au <- au3[, c("rigged1", "rigged2", "rigged3", "rigged4", "rigged5", "rigged6")]
item_intercor(interitem_au)
interitem_au <- ch3[, c("rigged1", "rigged2", "rigged3", "rigged4", "rigged5", "rigged6")]
item_intercor(interitem_au)
interitem_au <- de3[, c("rigged1", "rigged2", "rigged3", "rigged4", "rigged5", "rigged6")]
item_intercor(interitem_au)
interitem_au <- fr3[, c("rigged1", "rigged2", "rigged3", "rigged4", "rigged5", "rigged6")]
item_intercor(interitem_au)
interitem_au <- uk3[, c("rigged1", "rigged2", "rigged3", "rigged4", "rigged5", "rigged6")]
item_intercor(interitem_au)
interitem_au <- us3[, c("rigged1", "rigged2", "rigged3", "rigged4", "rigged5", "rigged6")]
item_intercor(interitem_au)



###################################
###################################
####### Online Appendix E  ########
###################################
###################################

# Differences in Means

diff_au_rigged # Australia
diff_fr_rigged # France
diff_de_rigged # Germany
diff_ch_rigged # Switzerland
diff_uk_rigged # UK
diff_us_rigged # US

# Lower CIs

diff_au_rigged - (1.64 * diff_au_rigged_se) # Australia
diff_fr_rigged - (1.64 * diff_fr_rigged_se) # France
diff_de_rigged - (1.64 * diff_de_rigged_se) # Germany
diff_ch_rigged - (1.64 * diff_ch_rigged_se) # Switzerland
diff_uk_rigged - (1.64 * diff_uk_rigged_se) # UK
diff_us_rigged - (1.64 * diff_us_rigged_se) # US

# Upper CIs

diff_au_rigged + (1.64 * diff_au_rigged_se) # Australia
diff_fr_rigged + (1.64 * diff_fr_rigged_se) # France
diff_de_rigged + (1.64 * diff_de_rigged_se) # Germany
diff_ch_rigged + (1.64 * diff_ch_rigged_se) # Switzerland
diff_uk_rigged + (1.64 * diff_uk_rigged_se) # UK
diff_us_rigged + (1.64 * diff_us_rigged_se) # US




###################################
###################################
####### Online Appendix F  ########
###################################
###################################

##### Figure F1: Treatment Effect on Redistributive Preference by Country Using Unweighted Data (applying the weighting variable NOWeights will cause the program to assume weights for y are uniformly 1)


diff_au_rigged <- -wtd.t.test(au_nb$redist_avg,au_rig$redist_avg,weight=au_nb$NOWEIGHTS,weighty=au_rig$NOWEIGHTS,samedata=F)$additional[1]
diff_au_rigged_se <- wtd.t.test(au_nb$redist_avg,au_rig$redist_avg,weight=au_nb$NOWEIGHTS,weighty=au_rig$NOWEIGHTS,samedata=F)$additional[4]

diff_fr_rigged <- -wtd.t.test(fr_nb$redist_avg,fr_rig$redist_avg,weight=fr_nb$NOWEIGHTS,weighty=fr_rig$NOWEIGHTS,samedata=F)$additional[1]
diff_fr_rigged_se <- wtd.t.test(fr_nb$redist_avg,fr_rig$redist_avg,weight=fr_nb$NOWEIGHTS,weighty=fr_rig$NOWEIGHTS,samedata=F)$additional[4]

diff_de_rigged <- -wtd.t.test(de_nb$redist_avg,de_rig$redist_avg,weight=de_nb$NOWEIGHTS,weighty=de_rig$NOWEIGHTS,samedata=F)$additional[1]
diff_de_rigged_se <- wtd.t.test(de_nb$redist_avg,de_rig$redist_avg,weight=de_nb$NOWEIGHTS,weighty=de_rig$NOWEIGHTS,samedata=F)$additional[4]

diff_ch_rigged <- -wtd.t.test(ch_nb$redist_avg,ch_rig$redist_avg,weight=ch_nb$NOWEIGHTS,weighty=ch_rig$NOWEIGHTS,samedata=F)$additional[1]
diff_ch_rigged_se <- wtd.t.test(ch_nb$redist_avg,ch_rig$redist_avg,weight=ch_nb$NOWEIGHTS,weighty=ch_rig$NOWEIGHTS,samedata=F)$additional[4]

diff_uk_rigged <- -wtd.t.test(uk_nb$redist_avg,uk_rig$redist_avg,weight=uk_nb$NOWEIGHTS,weighty=uk_rig$NOWEIGHTS,samedata=F)$additional[1]
diff_uk_rigged_se <- wtd.t.test(uk_nb$redist_avg,uk_rig$redist_avg,weight=uk_nb$NOWEIGHTS,weighty=uk_rig$NOWEIGHTS,samedata=F)$additional[4]

diff_us_rigged <- -wtd.t.test(us_nb$redist_avg,us_rig$redist_avg,weight=us_nb$NOWEIGHTS,weighty=us_rig$NOWEIGHTS,samedata=F)$additional[1]
diff_us_rigged_se <- wtd.t.test(us_nb$redist_avg,us_rig$redist_avg,weight=us_nb$NOWEIGHTS,weighty=us_rig$NOWEIGHTS,samedata=F)$additional[4]


# Plot

par(mar=c(6.1,5.1,2,2.1))
plot(x=c(),y=c(),ylim=c(-0.3,0.3),xlim=c(0.1,0.8), ylab="Treatment Effect", main="", xaxt="n", xlab="",cex.lab=1.4)
points(0.2,diff_au_rigged,pch=16)
lines(x=c(0.2,0.2),y=c(diff_au_rigged-1.64*diff_au_rigged_se,diff_au_rigged+1.64*diff_au_rigged_se),lwd=2)
points(0.3,diff_fr_rigged,pch=16)
lines(x=c(0.3,0.3),y=c(diff_fr_rigged-1.64*diff_fr_rigged_se,diff_fr_rigged+1.64*diff_fr_rigged_se),lwd=2)
points(0.4,diff_de_rigged,pch=16)
lines(x=c(0.4,0.4),y=c(diff_de_rigged-1.64*diff_de_rigged_se,diff_de_rigged+1.64*diff_de_rigged_se),lwd=2)
points(0.5,diff_ch_rigged,pch=16)
lines(x=c(0.5,0.5),y=c(diff_ch_rigged-1.64*diff_ch_rigged_se,diff_ch_rigged+1.64*diff_ch_rigged_se),lwd=2)
points(0.6,diff_uk_rigged,pch=16)
lines(x=c(0.6,0.6),y=c(diff_uk_rigged-1.64*diff_uk_rigged_se,diff_uk_rigged+1.64*diff_uk_rigged_se),lwd=2)
points(0.7,diff_us_rigged,pch=16)
lines(x=c(0.7,0.7),y=c(diff_us_rigged-1.64*diff_us_rigged_se,diff_us_rigged+1.64*diff_us_rigged_se),lwd=2)
abline(h=0,lty="dashed")
axis(1, at=seq(0.2,0.7,by=0.1), labels = FALSE)
text(seq(0.03,0.73,by=0.1), par("usr")[3] - 0.10, labels = c("","Australia","France","Germany","Switzerland", "United Kingdom","United States",""), srt = 45, pos = 1, xpd = TRUE,cex=1.2)



##### Figure F2: Treatment Effect on Redistributive Preference by Country with the Dependent Variable Constructed Using Factor Scores

# Save numbers from the models:

diff_au_rigged <- -wtd.t.test(au_nb$redist_fac_1to5,au_rig$redist_fac_1to5,weight=au_nb$weight3,weighty=au_rig$weight3,samedata=F)$additional[1]
diff_au_rigged_se <- wtd.t.test(au_nb$redist_fac_1to5,au_rig$redist_fac_1to5,weight=au_nb$weight3,weighty=au_rig$weight3,samedata=F)$additional[4]

diff_fr_rigged <- -wtd.t.test(fr_nb$redist_fac_1to5,fr_rig$redist_fac_1to5,weight=fr_nb$weight3,weighty=fr_rig$weight3,samedata=F)$additional[1]
diff_fr_rigged_se <- wtd.t.test(fr_nb$redist_fac_1to5,fr_rig$redist_fac_1to5,weight=fr_nb$weight3,weighty=fr_rig$weight3,samedata=F)$additional[4]

diff_de_rigged <- -wtd.t.test(de_nb$redist_fac_1to5,de_rig$redist_fac_1to5,weight=de_nb$weight3,weighty=de_rig$weight3,samedata=F)$additional[1]
diff_de_rigged_se <- wtd.t.test(de_nb$redist_fac_1to5,de_rig$redist_fac_1to5,weight=de_nb$weight3,weighty=de_rig$weight3,samedata=F)$additional[4]

diff_ch_rigged <- -wtd.t.test(ch_nb$redist_fac_1to5,ch_rig$redist_fac_1to5,weight=ch_nb$weight3,weighty=ch_rig$weight3,samedata=F)$additional[1]
diff_ch_rigged_se <- wtd.t.test(ch_nb$redist_fac_1to5,ch_rig$redist_fac_1to5,weight=ch_nb$weight3,weighty=ch_rig$weight3,samedata=F)$additional[4]

diff_uk_rigged <- -wtd.t.test(uk_nb$redist_fac_1to5,uk_rig$redist_fac_1to5,weight=uk_nb$weight3,weighty=uk_rig$weight3,samedata=F)$additional[1]
diff_uk_rigged_se <- wtd.t.test(uk_nb$redist_fac_1to5,uk_rig$redist_fac_1to5,weight=uk_nb$weight3,weighty=uk_rig$weight3,samedata=F)$additional[4]

diff_us_rigged <- -wtd.t.test(us_nb$redist_fac_1to5,us_rig$redist_fac_1to5,weight=us_nb$weight3,weighty=us_rig$weight3,samedata=F)$additional[1]
diff_us_rigged_se <- wtd.t.test(us_nb$redist_fac_1to5,us_rig$redist_fac_1to5,weight=us_nb$weight3,weighty=us_rig$weight3,samedata=F)$additional[4]


# Plot

par(mar=c(6.1,5.1,2,2.1))
plot(x=c(),y=c(),ylim=c(-0.3,0.3),xlim=c(0.1,0.8), ylab="Treatment Effect", main="", xaxt="n", xlab="",cex.lab=1.4)
points(0.2,diff_au_rigged,pch=16)
lines(x=c(0.2,0.2),y=c(diff_au_rigged-1.64*diff_au_rigged_se,diff_au_rigged+1.64*diff_au_rigged_se),lwd=2)
points(0.3,diff_fr_rigged,pch=16)
lines(x=c(0.3,0.3),y=c(diff_fr_rigged-1.64*diff_fr_rigged_se,diff_fr_rigged+1.64*diff_fr_rigged_se),lwd=2)
points(0.4,diff_de_rigged,pch=16)
lines(x=c(0.4,0.4),y=c(diff_de_rigged-1.64*diff_de_rigged_se,diff_de_rigged+1.64*diff_de_rigged_se),lwd=2)
points(0.5,diff_ch_rigged,pch=16)
lines(x=c(0.5,0.5),y=c(diff_ch_rigged-1.64*diff_ch_rigged_se,diff_ch_rigged+1.64*diff_ch_rigged_se),lwd=2)
points(0.6,diff_uk_rigged,pch=16)
lines(x=c(0.6,0.6),y=c(diff_uk_rigged-1.64*diff_uk_rigged_se,diff_uk_rigged+1.64*diff_uk_rigged_se),lwd=2)
points(0.7,diff_us_rigged,pch=16)
lines(x=c(0.7,0.7),y=c(diff_us_rigged-1.64*diff_us_rigged_se,diff_us_rigged+1.64*diff_us_rigged_se),lwd=2)
abline(h=0,lty="dashed")
axis(1, at=seq(0.2,0.7,by=0.1), labels = FALSE)
text(seq(0.03,0.73,by=0.1), par("usr")[3] - 0.10, labels = c("","Australia","France","Germany","Switzerland", "United Kingdom","United States",""), srt = 45, pos = 1, xpd = TRUE,cex=1.2)





###################################
###################################
####### Online Appendix G  ########
###################################
###################################


######## Five Country Bloc ########

model.m <- lm(elitedom_polychoric_factor_score ~ treatment + age + female + income + educ_fourlevels + lr + country, data = subset(combo, FiveCountries == "Yes"))
model.y <- lm(redist_avg ~ elitedom_polychoric_factor_score + treatment + age + female + income + educ_fourlevels + lr + country, data = subset(combo, FiveCountries == "Yes"))

out.1 <- mediate(model.m, model.y, sims = 1000, boot = TRUE, treat = "treatment", mediator = "elitedom_polychoric_factor_score")


summary(model.m) # Result of Model 1 - M1 
summary(model.y) # Result of Model 1 - Y1 
summary(out.1)   # ACME, direct effect, total effect


######## United States ########

model.m <- lm(elitedom_polychoric_factor_score ~ treatment + age + female + income + educ_fourlevels + lr, data = subset(combo, country == "US"))
model.y <- lm(redist_avg ~ elitedom_polychoric_factor_score + treatment + age + female + income + educ_fourlevels + lr, data = subset(combo, country == "US"))

out.1 <- mediate(model.m, model.y, sims = 1000, boot = TRUE, treat = "treatment", mediator = "elitedom_polychoric_factor_score")

summary(model.m) # Result of Model 1 - M1 
summary(model.y) # Result of Model 1 - Y1 
summary(out.1) # ACME, direct effect, total effect




###################################
###################################
####### Online Appendix H  ########
###################################
###################################



############ Step 1 – First stage results ############


# United States, with covariates
Step.One.Mediation <- lm(elitedom_polychoric_factor_score ~ treatment 
                         + age + female + as.numeric(educ_fourlevels) + income + lr, data = subset(combo, country == "US"))
summary(Step.One.Mediation)

# Five-Country Boc, with covariates
Step.One.Mediation <- lm(elitedom_polychoric_factor_score ~ treatment 
                         + age + female + as.numeric(educ_fourlevels) + income + lr, data = subset(combo, country != "US"))
summary(Step.One.Mediation)


# United States, without covariates
Step.One.Mediation <- lm(elitedom_polychoric_factor_score ~ treatment 
                         , data = subset(combo, country == "US"))
summary(Step.One.Mediation)

# Five-Country Boc, without covariates
Step.One.Mediation <- lm(elitedom_polychoric_factor_score ~ treatment 
                         , data = subset(combo, country != "US"))
summary(Step.One.Mediation)


############ Step 2 – Reduced form results ############

# United States, with covariates
Step.Two.Mediation <- lm(redist_avg ~ treatment 
                         + age + female + as.numeric(educ_fourlevels) + income + lr, data = subset(combo, country == "US"))
summary(Step.Two.Mediation)

# Five-Country Boc, with covariates
Step.Two.Mediation <- lm(redist_avg ~ treatment 
                         + age + female + as.numeric(educ_fourlevels) + income + lr, data = subset(combo, country != "US"))
summary(Step.Two.Mediation)

# United States, without covariates
Step.Two.Mediation <- lm(redist_avg ~ treatment 
                         , data = subset(combo, country == "US"))
summary(Step.Two.Mediation)

# Five-Country Boc, without covariates
Step.Two.Mediation <- lm(redist_avg ~ treatment 
                         , data = subset(combo, country != "US"))
summary(Step.Two.Mediation)


############ Step 3 – 2SLS results ############

# United States, with covariates
Step.Three.Mediation <- ivreg(redist_avg ~ elitedom_polychoric_factor_score + age + female + as.numeric(educ_fourlevels) + income + lr
                              | treatment + age + female + as.numeric(educ_fourlevels) + income + lr,  data = subset(combo, country == "US"))
summary(Step.Three.Mediation)

# Five-Country Boc, with covariates
Step.Three.Mediation <- ivreg(redist_avg ~ elitedom_polychoric_factor_score + age + female + as.numeric(educ_fourlevels) + income + lr
                              | treatment + age + female + as.numeric(educ_fourlevels) + income + lr,  data = subset(combo, country != "US"))
summary(Step.Three.Mediation)


# United States, without covariates
Step.Three.Mediation <- ivreg(redist_avg ~ elitedom_polychoric_factor_score 
                              | treatment ,  data = subset(combo, country == "US"))
summary(Step.Three.Mediation)

# Five-Country Boc, without covariates
Step.Three.Mediation <- ivreg(redist_avg ~ elitedom_polychoric_factor_score 
                              | treatment ,  data = subset(combo, country != "US"))
summary(Step.Three.Mediation)





###################################
###################################
####### Online Appendix I  ########
###################################
###################################

summary(PredictingAU_stand)
summary(PredictingCH_stand)
summary(PredictingDE_stand)
summary(PredictingFR_stand)
summary(PredictingUK_stand)
summary(PredictingUS_stand)

# To print out the regression tables in a combined format

stargazer(PredictingAU_stand, PredictingCH_stand, PredictingDE_stand,
          title="columnsAlternative", align=TRUE,
          type = "html",
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "Table_I1_part1.doc")

stargazer(PredictingFR_stand, PredictingUK_stand, PredictingUS_stand,
          title="columnsAlternative", align=TRUE,
          type = "html",
          star.cutoffs = c(0.05, 0.01, 0.001),
          out = "Table_I1_part2.doc")



###################################
###################################
####### Online Appendix J ########
###################################
###################################

# Model 5a
mod2 <- lm(redist_avg ~ treatment*lr, data=subset(combo, country=="US"), weights = weight3)
summary(mod2)

# Model 5b
mod2 <- lm(redist_avg ~ treatment*lr, data=subset(combo, FiveCountries=="Yes"), weights = weight3)
summary(mod2)

# Model K1a
mod2 <- lm(redist_avg ~ treatment*thermo_diff, data=subset(combo, country=="US"), weights = weight3)
summary(mod2)

# Model K1b
mod2 <- lm(redist_avg ~ treatment*thermo_diff, data=subset(combo, FiveCountries=="Yes"), weights = weight3)
summary(mod2)

# Model 7a
mod2 <- lm(redist_avg ~ treatment*lf4, data=subset(combo, country=="US"), weights = weight3)
summary(mod2)

# Model 7b
mod2 <- lm(redist_avg ~ treatment*lf4, data=subset(combo, FiveCountries=="Yes"), weights = weight3)
summary(mod2)

# Model K2a
mod2 <- lm(redist_avg ~ treatment*rr_avg, data=subset(combo, country=="US"), weights = weight3)
summary(mod2)

# Model K2b
mod2 <- lm(redist_avg ~ treatment*rr_avg, data=subset(combo, FiveCountries=="Yes"), weights = weight3)
summary(mod2)

# Model K4a
mod2 <- lm(redist_avg ~ treatment*indiv_avg, data=subset(combo, country=="US"), weights = weight3)
summary(mod2)

# Model K4b
mod2 <- lm(redist_avg ~ treatment*indiv_avg, data=subset(combo, FiveCountries=="Yes"), weights = weight3)
summary(mod2)






###################################
###################################
####### Online Appendix K ########
###################################
###################################

# Figure K1a

mod2 <- lm(redist_avg ~ treatment*thermo_diff, data=subset(combo, country=="US"), weights = weight3)
summary(mod2)

par(mfrow=c(1,1),oma=c(1,1,0,0),mai=c(0.5,0.4,0,0),mar=c(4,4,1,1)) # oma is originally 0 0 0 0. mai is 0.6732 0.5412 0.5412 0.2772
thermo_sim <- -10:10

treat <- coef(mod2)[2] + coef(mod2)[4] * thermo_sim
treat_se <- sqrt(vcov(mod2)[2, 2] + thermo_sim^2 * vcov(mod2)[4, 4] + 2 * thermo_sim * vcov(mod2)[2, 4])
plot(x = thermo_sim, y = treat, type = "p",pch=16,
     xlab = "Partisan Identification",
     ylab = "Treatment Effect",
     ylim = c(-0.4,0.4), cex.lab=1.1)
for(i in 1:21){
  
  lines(c(i-11,i-11),c(treat[i],treat[i] + 1.96 * treat_se[i]),lwd=2)
  lines(c(i-11,i-11),c(treat[i],treat[i] - 1.96 * treat_se[i]),lwd=2)
  
}
abline(h=0,lty="dashed")

par(new=T)
hist(subset(combo, country == "US")$thermo_diff,xlab=NULL,ylab=NULL,axes=F,main=NULL,freq=F,col=rgb(0.1,0.1,0.1,0.1),ylim=c(0,1.0))


# Figure K1b

mod2 <- lm(redist_avg ~ treatment*thermo_diff, data=subset(combo, FiveCountries=="Yes"), weights = weight3)
summary(mod2)


par(mfrow=c(1,1),oma=c(1,1,0,0),mai=c(0.5,0.4,0,0),mar=c(4,4,1,1)) # oma is originally 0 0 0 0. mai is 0.6732 0.5412 0.5412 0.2772
thermo_sim <- -10:10

treat <- coef(mod2)[2] + coef(mod2)[4] * thermo_sim
treat_se <- sqrt(vcov(mod2)[2, 2] + thermo_sim^2 * vcov(mod2)[4, 4] + 2 * thermo_sim * vcov(mod2)[2, 4])
plot(x = thermo_sim, y = treat, type = "p",pch=16,
     xlab = "Partisan Identification",
     ylab = "Treatment Effect",
     ylim = c(-0.4,0.4), cex.lab=1.1)
for(i in 1:21){
  
  lines(c(i-11,i-11),c(treat[i],treat[i] + 1.96 * treat_se[i]),lwd=2)
  lines(c(i-11,i-11),c(treat[i],treat[i] - 1.96 * treat_se[i]),lwd=2)
  
}
abline(h=0,lty="dashed")

par(new=T)
hist(subset(combo, FiveCountries == "Yes")$thermo_diff,xlab=NULL,ylab=NULL,axes=F,main=NULL,freq=F,col=rgb(0.1,0.1,0.1,0.1),ylim=c(0,1.0))



# Figure K2a

mod2 <- lm(redist_avg ~ treatment*rr_avg, data=subset(combo, country=="US"), weights = weight3)
summary(mod2)


par(mfrow=c(1,1),oma=c(1,1,0,0),mai=c(0.5,0.4,0,0),mar=c(4,4,1,1)) # oma is originally 0 0 0 0. mai is 0.6732 0.5412 0.5412 0.2772
thermo_sim <- seq(from = 1, to = 5, by = .25)

treat <- coef(mod2)[2] + coef(mod2)[4] * thermo_sim
treat_se <- sqrt(vcov(mod2)[2, 2] + thermo_sim^2 * vcov(mod2)[4, 4] + 2 * thermo_sim * vcov(mod2)[2, 4])
plot(x = thermo_sim, y = treat, type = "p",pch=16,
     xlab = "Racial Resentment",
     ylab = "Treatment Effect",
     ylim = c(-1.0,1.2), cex.lab=1.1)
for(i in 1:17){
  lines(c(((i/i)+(i*.25)-(.25)),((i/i)+(i*.25)-(.25))),c(treat[i],treat[i] + 1.96 * treat_se[i]),lwd=2)
  lines(c(((i/i)+(i*.25)-(.25)),((i/i)+(i*.25)-(.25))),c(treat[i],treat[i] - 1.96 * treat_se[i]),lwd=2)
}
abline(h=0,lty="dashed")

par(new=T)
hist(subset(combo, country == "US")$rr_avg,xlab=NULL,ylab=NULL,axes=F,main=NULL,freq=F,col=rgb(0.1,0.1,0.1,0.1),ylim=c(0,2.5))



# Figure K2b

mod2 <- lm(redist_avg ~ treatment*rr_avg, data=subset(combo, FiveCountries=="Yes"), weights = weight3)
summary(mod2)


par(mfrow=c(1,1),oma=c(1,1,0,0),mai=c(0.5,0.4,0,0),mar=c(4,4,1,1)) # oma is originally 0 0 0 0. mai is 0.6732 0.5412 0.5412 0.2772
thermo_sim <- seq(from = 1, to = 5, by = .25)

treat <- coef(mod2)[2] + coef(mod2)[4] * thermo_sim
treat_se <- sqrt(vcov(mod2)[2, 2] + thermo_sim^2 * vcov(mod2)[4, 4] + 2 * thermo_sim * vcov(mod2)[2, 4])
plot(x = thermo_sim, y = treat, type = "p",pch=16,
     xlab = "Racial Resentment",
     ylab = "Treatment Effect",
     ylim = c(-1.0,1.2), cex.lab=1.1)
for(i in 1:17){
  lines(c(((i/i)+(i*.25)-(.25)),((i/i)+(i*.25)-(.25))),c(treat[i],treat[i] + 1.96 * treat_se[i]),lwd=2)
  lines(c(((i/i)+(i*.25)-(.25)),((i/i)+(i*.25)-(.25))),c(treat[i],treat[i] - 1.96 * treat_se[i]),lwd=2)
}
abline(h=0,lty="dashed")

par(new=T)
hist(subset(combo, FiveCountries == "Yes")$rr_avg,xlab=NULL,ylab=NULL,axes=F,main=NULL,freq=F,col=rgb(0.1,0.1,0.1,0.1),ylim=c(0,2.5))




# Figure K3

diff_us_rigged_white_all <- -wtd.t.test(subset(us_nb, White == "White")$redist_avg,
                                        subset(us_rig, White == "White")$redist_avg,
                                        #                               weight=subset(us_nb, white == "White")$weight3,
                                        #                              weighty=subset(us_rig, white == "White")$weight3,
                                        samedata=F)$additional[1]


diff_us_rigged_se_white_all <- wtd.t.test(subset(us_nb, White == "White")$redist_avg,
                                          subset(us_rig, White == "White")$redist_avg,
                                          #                                   weight=subset(us_nb, white == "White")$weight3,
                                          #                                  weighty=subset(us_rig, white == "White")$weight3,
                                          samedata=F)$additional[4]

diff_us_rigged_nonwhite_all <- -wtd.t.test(subset(us_nb, White == "Minority")$redist_avg,
                                           subset(us_rig, White == "Minority")$redist_avg,
                                           #                               weight=subset(us_nb, white == "Minority")$weight3,
                                           #                              weighty=subset(us_rig, white == "Minority")$weight3,
                                           samedata=F)$additional[1]


diff_us_rigged_se_nonwhite_all <- wtd.t.test(subset(us_nb, White == "Minority")$redist_avg,
                                             subset(us_rig, White == "Minority")$redist_avg,
                                             #                                weight=subset(us_nb, white == "Minority")$weight3,
                                             #                               weighty=subset(us_rig, white == "Minority")$weight3,
                                             samedata=F)$additional[4]



#pdf("plot_w3_effects_regulfac1to5.pdf",width=8)
par(mar=c(3.1,5.1,2,2.1))
plot(x=c(),y=c(),ylim=c(-0.5,0.5),xlim=c(0.1,0.4), ylab="Treatment Effect", main="", xaxt="n", xlab="",cex.lab=1.4)

points(0.2,diff_us_rigged_white_all,pch=16)
lines(x=c(0.2,0.2),y=c(diff_us_rigged_white_all-1.96*diff_us_rigged_se_white_all,diff_us_rigged_white_all+1.96*diff_us_rigged_se_white_all),lwd=2)
points(0.3,diff_us_rigged_nonwhite_all,pch=16, col = "black")
lines(x=c(0.3,0.3),y=c(diff_us_rigged_nonwhite_all-1.96*diff_us_rigged_se_nonwhite_all,diff_us_rigged_nonwhite_all+1.96*diff_us_rigged_se_nonwhite_all),lwd=2, col = "black")

abline(h=0,lty="dashed")
axis(1, at=seq(0.2,0.3,by=0.1), labels = FALSE)
text(seq(0.1,0.8,by=0.1), par("usr")[3] - 0.03, labels = c("","White","Minority"), srt = 0, pos = 1, xpd = TRUE,cex=1.0)
#dev.off()



# Figure K4a

mod2 <- lm(redist_avg ~ treatment*indiv_avg, data=subset(combo, country=="US"), weights = weight3)
summary(mod2)

par(mfrow=c(1,1),oma=c(1,1,0,0),mai=c(0.5,0.4,0,0),mar=c(4,4,3,1)) # oma is originally 0 0 0 0. mai is 0.6732 0.5412 0.5412 0.2772
#thermo_sim <- 1:5
thermo_sim <- seq(from = 1, to = 5, by = .20)


treat <- coef(mod2)[2] + coef(mod2)[4] * thermo_sim
treat_se <- sqrt(vcov(mod2)[2, 2] + thermo_sim^2 * vcov(mod2)[4, 4] + 2 * thermo_sim * vcov(mod2)[2, 4])
plot(x = thermo_sim, y = treat, type = "p",pch=16,
     #  main = "Interaction of Treatment * Individualism - US",
     xlab = "Individualist Views",
     ylab = "Treatment Effect",
     ylim = c(-1.0,1.2), cex.lab=1.1)
for(i in 1:21){
  lines(c(((i/i)+(i*.2)-(.2)),((i/i)+(i*.2)-(.2))),c(treat[i],treat[i] + 1.96 * treat_se[i]),lwd=2)
  lines(c(((i/i)+(i*.2)-(.2)),((i/i)+(i*.2)-(.2))),c(treat[i],treat[i] - 1.96 * treat_se[i]),lwd=2)
}
abline(h=0,lty="dashed")

par(new=T)
hist(subset(combo, country == "US")$indiv_avg,xlab=NULL,ylab=NULL,axes=F,main=NULL,freq=F,col=rgb(0.1,0.1,0.1,0.1),ylim=c(0,2.5), breaks = 21)


# Figure  K4b

mod2 <- lm(redist_avg ~ treatment*indiv_avg, data=subset(combo, FiveCountries=="Yes"), weights = weight3)
summary(mod2)


par(mfrow=c(1,1),oma=c(1,1,0,0),mai=c(0.5,0.4,0,0),mar=c(4,4,3,1)) # oma is originally 0 0 0 0. mai is 0.6732 0.5412 0.5412 0.2772
#thermo_sim <- 1:5
thermo_sim <- seq(from = 1, to = 5, by = .20)


treat <- coef(mod2)[2] + coef(mod2)[4] * thermo_sim
treat_se <- sqrt(vcov(mod2)[2, 2] + thermo_sim^2 * vcov(mod2)[4, 4] + 2 * thermo_sim * vcov(mod2)[2, 4])
plot(x = thermo_sim, y = treat, type = "p",pch=16,
     #  main = "Interaction of Treatment * Individualism - US",
     xlab = "Individualist Views",
     ylab = "Treatment Effect",
     ylim = c(-1.0,1.2), cex.lab=1.1)
for(i in 1:21){
  lines(c(((i/i)+(i*.2)-(.2)),((i/i)+(i*.2)-(.2))),c(treat[i],treat[i] + 1.96 * treat_se[i]),lwd=2)
  lines(c(((i/i)+(i*.2)-(.2)),((i/i)+(i*.2)-(.2))),c(treat[i],treat[i] - 1.96 * treat_se[i]),lwd=2)
}
abline(h=0,lty="dashed")

par(new=T)
hist(subset(combo, FiveCountries == "Yes")$indiv_avg,xlab=NULL,ylab=NULL,axes=F,main=NULL,freq=F,col=rgb(0.1,0.1,0.1,0.1),ylim=c(0,2.5), breaks = 20)


